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ABSTRACT 


The three-dimensional, primitive equations of motion have been solved 
numerically for the case of isotropic box turbulence and the distortion of 
homogeneous turbulence by irrotational plane strain at large Reynolds 
numbers. A Gaussian filter was applied to governing equations to define 
the large scale field. This gives rise to additional second order computed 
scale stresses (Leonard stresses). The residual stresses are simulated 
through an eddy viscosity. 16x16x16 and 32x32x32 uniform grids were used, 
with a fourth order differencing scheme in space and a second order 
Adams-Bashforth predictor for explicit time stepping. The results were 
compared to the experiments and statistical information was extracted 
from the computer generated data. 
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CHAPTER I 


INTRODUCTION 


1.1 Historical Background 

As computer capabilities grow, the three-dimensional time-dependent 
computation of turbulence is becoming possible. However, the retention 
of all scales of motion is not yet feasible (and probably never will be), 
so the best one can hope for is the simulation of the large scale 
structures. The large scale structures are strongly dependent upon the 
nature of the flow, but there is considerable evidence that the structure 
of the smaller scales is independent of the large scale structure. This 
suggests a mixed approach in which one computes the large scale motions 
and models the small scales . 

To define the large scale motions, some sort of averaging operator 
has to be applied to the governing equations to filter out the small 
scale motions. However, in three-dimensional computations to date, the 
definition of the small scale motions was not precisely related to the 
filtering operation and consequently the meaning of those motions was 
not very clear. In any case, the resulting equations for the filtered 
field contain so-called sub-grid scale or residual scale Reynolds stresses, 
which must be modeled in the computation. 

Two distinctive solution methods have been used in solving the re- 
sulting equations. The first is a conventional finite difference mesh 
calculation. In this approach, the simplest and perhaps the most usual 
way of relating the residual scale turbulence to the filtered motion is 
by a local eddy viscosity model. Smagorinsky (1963, 1965) related eddy 
viscosity to the local strain-rate of the filtered field. Deardorff 
(1970a, b) applied this model to three-dimensional turbulent channel flow 
and planetary boundary layer problems, and Schumann (1973) applied it to 
plane channels and annuli. Deardorff (1973) and Schumann later introduced 
more sophisticated residual Reynolds stress transport equations. Other 
examples of this approach are given by Lilly (1964, 1967), Smagorinsky 
et al (1965) , Fox and Lilly (1972) , Fox and Deardorff (1972) . Previous 
work mentioned above has not paid sufficient attention to the basic 
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aspects of this type of simulation, so as yet this approach has not 
really progressed very far. 

The second approach is the spectral (Fourier) method much advocated 
by Orszag (1969, 1971a, b) . While this method has mathematically attrac- 
tive features for certain problems, it is generally more difficult to 
extend to flows with interesting geometries. Moreover, work to date 
ignored the residual Reynolds stresses, and it is not clear how these 
could be incorporated in a Fourier calculation. 

1.2 Motivation and Objectives 

At the time this work was initiated there were some serious problems 
with work done previously: 

(1) There was a need to define the large-scale field precisely, 
so that the equations can be systematically developed. 

(2) There was a need to carefully evaluate the accuracy requirements 
to be sure that computational errors are higher order than 

the residual stresses. 

(3) There was a need to carefully assess just what can really 
be learned from this type of turbulence simulation. 

The main objective of this study was to carefully develop a numerical 
simulation method for turbulent flows away from solid or free boundaries, 
and to apply this method to study decaying isotropic turbulence and 
homogeneous turbulence with irrotational plane strain. 

The present study is one in a systematic program investigating large 
eddy simulation of turbulence, and reports the details of the initial 
computations under this program. 

1.3 Summary 

The contribution of the present work includes 

(a) a precise definition for the large-scale field (after 
Leonard (1973)), 

(b) a study of the optimum averaging scale, as compared to the 
grid mesh scale, 

(c) a study of two residual stress models, and evaluation of the 
model constant for each, 
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(d) a demonstration that the model constant Is independent of 
mesh size, 

(e) a fourth-order differencing scheme that properly conserves 
energy and momentum, 

(f) a method for calculating the pressure so as to conserve 
mass at subsequent time steps, 

(g) a demonstration that a coarse mesh can be used to obtain 
surprisingly good predictions for the Reynolds stresses in 
a straining flow, 

(h) an evaluation of certain aspects of simple turbulence closure 
models . 

Although many questions have been answered by this work, new ones 
have been raised. Suggestions for follow-on work in this project are 
made in Chapter VI. 
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Chapter II 


Theoretical Foundations 

2 . 1 Definition of the Filtered and Residual Fields 

To resolve the smallest scales of turbulence in a grid-based cal- 
culation, the mesh size has to be smaller than the dissipation length, 

3 1/4 

which is on the order of the Kolmogorov microscale, r| - (v /e) 

Here v is the kinematic viscosity and e is the energy dissipation 

rate per unit mass. It is known that (see Tennekes and Lumley 1972) 

3 

e - q /L where q is the R.M.S. velocity and L is the length scale 
of large eddies. Thus the minimum number of mesh points that must be 
used in a three-dimensional grid computation that resolves both the 
large and small scales can be estimated as 



3 

Using this estimate, we find that an R_ of 10 , typical of turbulent 

7 1 

flows, would require 2 x 10 worlds of storage for four variables. This 
is approximately 50 times the Large Core Memory of the CDC 7600, and 
about 10 times the available memory of the ILLIAC IV disk. 

It is clear that one can not do a full simulation, except at 
extremely low Reynolds number. The best one can hope for is a computa- 
tion that will yield the large-scale motions . Fortunately these contain 
most of the turbulence energy, and are responsible for most of the tur- 
bulent transport, and so a large-eddy simulation technique would be 
very useful, especially if it could handle arbitrary flows. 

The first problem one faces is in defining the large scale field. 
Conventionally the large scale motions have been defined by volume- 
averaging in a continuous manner over computational grid boxes (e.g. 
Deardorff 1970a). Schumann (1973, 1974) applied a slightly different 

technique involving averaging over the surface of grid boxes. 

A more general approach that recognizes the continuous nature of 

the flow variables is the "filter function" approach of Leonard (1973) . 
Let f(x) denote a field variable, for example velocity, f may contain 
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large and very small components. Then, we define the filtered field 
f by 

where G is a selected filter function. For C ■ C , where C Is a 
constant, the filter G must satisfy 

G(x) dsc = 1 (2.1b) 

Now f can be decomposed into its filtered field (FF) component, f , 
and residual-field (RF) components, f* , by 

f » f + f ' (2.2) 



f(x) 


J” G (x~2i 


') f(x') dx* 


(2.1a) 


Note that f is not the conventional mean used in the classic-turbulence 
literature. 

In the present work we treat only flows that are homogeneous , for 
which the integration in (2.1) extends over all space. Careful considera- 
tion will have to be given to the domain of integration when one desires 
to treat a flow near a wall. 

Now note that, if G is piecewise continuously differentiable and 
G(r) goes to zero as r -*■ 00 at least as fast as 1/r^ , 


3f 

3x 


00 


/ 

—00 


G(x-x’) 


3f 

3x/ 


dx ? 


00 00 



-£0 


f (*’) G( £-£'> 
f(x') dx' 


dx' 



(2.3a) 


Also , 


3f _ 3f 

at " at 


(2.3b) 
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However, 


fl + fl 


(2.3c) 


2.2 The Dynamical Equations 

Applying (2.1) to the Navier-Stokes equations, and using (2.2) and 
(2.3), one obtains (for incompressible flows) 


3x i 

fj|± . — 

at 3Xj u i u j 


. I 3L + vy 2 0 

p 3x. i 


The advection term is 


u i u j “ u i u j + u i u j + u i u j + u i u j 


where 


" u i u j + R ij 


l ij " U I U j + U I*j + Vj 


(2.4) 

(2.5) 


(2.6) 


is the residual field contribution to the advection term. -pR^ 
called the "residual stress." 

To localize the first term on the right in (2.6) we carry out a Taylor 
series expansion, 


u lUj ( V t) 


-/ 


GOCq-x) VjUj (x) dx 
du 


J j u i<V t) + 3^ ‘W 5 + i ( V\o H V 




3u 


+ 0(x-x Q ) j > j Uj (^.t) (x k -x kQ ) 

1 3 } 

+ 23^ (: \-\o y + 0(x - x o ) j G( ^o-^ d * 

(2.7) 
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For the above filtering of the dynamical equations to be useful, the 
integrals in (2.7) must exist. This requires that G(r) -*• 0 exponentially 
as r 00 . 


2.3 Filter Selection 


(a) Sub-Grid-Scale filter 

Let us first seek a filter that makes the scales of motion in the 
residual field smaller than the scales in the filtered field, in the 
Fourier sense. Let 


00 

f - f f (k) dk 


( 2 . 8 ) 


Then, 


OO 00 


f - J J f(k) G(x-x') e 1 -- dx' dk 


(2.9) 


-«0 —00 


We want to have f contain all scales larger than a cut-off scale. Thus 
we want +k £ 

f - J f (k) e 1 -*- dk (2.10) 

-k 


where k is the cut-off wave number. Hence, we can write 
c ’ 


OO 00 00 

J H(k) f (k) e 1 -*- dk => J f (k) J G(x-x 


,v ik*x' . , ,, 

) e dx' dk 


where 


H(k;k c ) 


( 0 if k, > k 
1 c 

1 otherwise 


for any i 


(2.11a) 


(2.11b) 


So we have an integral equation for G , 

OO 


H(k; k^) » j" G(x-x') 


( 2 . 12 ) 


The solution to (2.12) is 


G(x-x’) 


3 sin ir(x.-x’)/A 

n 1 1 A 


i=l 


1T(Xi-Xp 


(2.13) 
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where A. = ir/k is the averaging or filtering length scale. This is 
the proper filter if one wishes to have the residual field really be 
"sub-grid scale." A grid-based computation made using this filter 
would be equivalent to a Fourier computation. 

The second moment of G involves integrals like 



sin(TTx/A A ) 

TPC 


dx 


(2.14) 


This integral does not exist, and hence the expansion (2.7) could not 
be used. This filter is not suitable for a grid-based numerical method; 
hence one can not expect to really have the residual field be sub-grid 
scale. 


(b) Top-hat filter 

The filter used implicitly by many workers is the top-hat. 


G(x-x') 


1/A, 


for x-x' < 


l£-i I 


(2.15) 


Then the filtered velocity is 

u (x) » 




A/ 2 


3 J «(*♦£> d£ 

A A - A A/2 


(2.16) 


This is equivalent to volume averaging. The Fourier transform of (2.16) is 

A A/2 

i r r 

u(k) 


f L 

t, l 

A A - A A/2 


jj (x+^) e — — dxdj^ 


-A A /2 

+ A A/2 


U(k) e 1 - ■ £ d£ 


{ 3 8in( 
i=l 


sin(k i A A/ / 2 ) 


Aa/2 


J 


(2.17) 


Here ^(k) is the Fourier transform of u . 
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Equation (2.17) shows that the spectrum of the filtered field will 
contain components of all wave-numbers. Moreover, at the wave-numbers 
for which the coefficient of uOO in (2.17) is zero, the inverse trans- 
form will be singular. This makes it impossible to predict the actual 

A 

spectrum ti(k) from the filtered spectrum u(k) , and this very unde- 
sirable feature of the top-hat filter renders it useless if we want to 
compute spectral features with a grid-based method. However, the top- 
hat filter could be used if spectral results were not sought. 

Using (2.15) in (2.7), and carrying out the integration over x , 
a 2 

Vj Q^.t) “ (^.t) + 2^ V 2 ( u i u j)+ O(A^) (2.18) 


2 2 - - 

The second term on the right is called the Leonard term; -pA^ V (u^u. )/24 
is called the "Leonard stress." As will be shown later, 


°“a> 


and A. 


0(A) 


So both the residual stresses and the 


ij “'“A" “““■ “A 

Leonard stresses have to be included; moreover, the computational dif- 

2 

ference scheme must be accurate to 0(A ) to avoid introduction of 
numerical errors comparable with these stresses. 


(c) Gaussain filter 

A filter with much more desirable properties is 


G(x-x') 



-Y(x-x') 2 /A a 


(2.19) 


where Y 


is a constant. Then the filtered velocity is 


u(x) 



-Y(x-x') 2 /A 2 

dx ' 


( 2 . 20 ) 


The Fourier transform of this is 
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u(x) - (VTr - ) 3 J f e_i “ “ e Y ? Ma d i 

' _oo _00 

■ (V?t) 3 / -® ei ~ -i e ’ 5 Ma di 

■ (VF^) “<« exp j- f- Y (^ + k 2 + k 3 > | 

" iiQc) ex P f - 4 y (2.21) 

Consider now the three-dimensional energy spectra of the actual and 
filtered fields: 


E(k) - 2 tt k 2 < G(k) • G*(k)> (2.22a) 

E(k) - 2 tt k 2 < u(k) • u* ( k ) > (2.22b) 

Here < > denotes an average over an ensemble of experiments, and 

* denotes a complex conjugate. Equations (2.21) and (2.22) show that 

( A l 2 ) 

E(k) = E(k) exp 1 - jj- k I (2.23) 

We see that the use of the Gaussian filter will result in a filtered 
field that misses only a very small amount of large scale motion; most 
of the small scale motions are placed in the residual field. Thus, in 
many respects this filter has the desirable properties of the sub-grid- 
scale filter. However, its behavior at r **■ 00 makes the integrals in 
(2.7) exist. Moreover, the conversion back and forth between the 
spectrum of the filtered field and the spectrum of the actual field is 

easily accomplished, and hence the Gaussian filter is perferable to the 

top-hat filter. 
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Using (2.19) and (2.7), one obtains 


+ 0(A*) (2.24) 

When y «* 6 , the Leonard term in (2.24) is exactly the same as in (2.18). 
Hence the Gaussian filter with y = 6 was chosen for the present study. 
This filter is illustrated on Fig. 2.1 and an example of E and E 
relation (2.23) is shown on Fig. 2.2. 

2.4 Residual Stress Models 

The following eddy viscosity model is used for 

R ij ■“ 3 \k 6 ij " 2 v t\j 



is the strain-rate tensor, and is an effective viscosity associated 

with the residual field motions. 

(a) Smagorinsky model 

Smagorinsky (1963) suggested a model for v T , 

V T - <c s A a ) 2 (2i y I^) 1 ' 2 (2.26) 

where Cg is a constant. In experiments one observes a sharp separation 
of turbulent regions, containing vorticity and non-turbulent regions which 
are irrotational. A weakness of this model is that, in a non-turbulent 
irrotational region, v T will have a non-zero value. This will give 
rise to residual stresses in the non-turbulent flow outside of a boundary 
layer . 


ij 


(2.25) 


u i u j ( v° 


“iVV 0 + 47 72< Vj > 
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(b) Vorticity model 

A way around this drawback is to relate directly to vorticity. 

A likely possibility is 

V T " (c v A A )2 V^X (2 ‘ 27 > 

where w - curl u is the vorticity, and c^ is a constant. 


2.5 Governing Equations for the Filtered Field 

Now, neglicting the molecular viscosity term, and dropping terms of 
2 

higher order than A^ , filtered momentum equations become 


3u 

It 


r + (vj + 2l v2 Vj - 2 Vij ) * - ^ 


where P ■ ^ + ^ R.. . This may be written as 
p 3 ii 


3 t 


. 3P A u 
h i " 3x^ " H i 


where 


l i " 3xj ( U i U j + 2 A 7 u i u j ” 2V T S ij) 


(2.28) 


(2.29) 


It is in this form that we shall deal with the problem computationally. 
The manner in which continuity was used to fix P is discussed in the 
next chapter. 
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CHAPTER III 


NUMERICAL METHOD 


3.1 Grid Layout and Notation 

A uniform cubic mesh is used, as sketched below. The mesh width A 
need not be the same as the averaging scale A^ introduced in the 
previous chapter. 



SL+ 1 


The i-component of the filtered velocity at the nth time step is 
written as 


(k,Z,m) 


where (k,£,m) is the meshpoint index for (x, y, z) . 

We now define the following operator notations: 

S/Sx.^ “ finite difference operator corresponding to 3/3x^ 

6/6t ■ finite difference operator corresponding to 3/3t 


G » finite difference form of gradient operator 
D ° finite difference form of divergence operator 


6 3 

(u^f) = transport operator corresponding to (u^f) 

Further details of these terms are given next . 
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3.2 Space Differencing 

A fourth order differencing scheme is applied where fourth order 
accuracy is needed. Since the Leonard and residual stress terms are 
second order, they can be approximated by second order formula to give 
the same accuracy. The central difference fourth order scheme is, for 
example, 


6u 

6x 


1 ) " 

12E | U (k-2) 


8u (k-l) + 8u (k+l) 


(k+2) 


(3.1) 


For simplicity, the subscripts £. and m are not shown. 

Suppose we represent u by a discrete Fourier expansion, 

u - T. u(k) e 1 -- (3.2) 

where 


n 


, 2ir 

k l ■ NA n l 


wave number in the 


*1 


direction 


n^ = -N/2 0,1,..., Oj - 1) 


N 


Number of mesh points in one direction 


The sum extends over all n^ , , and n^ . Substituting (3.2) in (3.1), 

the Fourier transform of 6u/6x is identified as 

. -i2Ak. -iAk.. iAk. i2Ak 
. 1 (e ■ 1 - 8e 1 + 8e 1 - . 1 ) u 


^ | 8 sln(Ak^) - Bin (2£k^ ) | u 


(3.3) 


If a modified wave number, k^ , is defined by 

k| = ^ | 8 sin(Ak^) - sin(2Ak^) ^ (3.4) 

then the Fourier transform of Du = 0 can be written as 

kp^ = 0 (3.5) 

A 

Note that the exact transform of div ij is k^u^ . Hence, k^ may be 
interpreted as the wave number that allows continuity to be satisfied in 
grid space. 
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If instead one were to use a second-order central difference scheme. 


6_u _1_ /- \ 

6x = 2A y (k+1) ” u (k-l ) j 


(3.6) 


The modified wave number , k^' , would be 


k^ = ^ sin(Ak i ) 


(3.7) 


k^, k^ and k^ are compared in Fig. 3.1. 

The fourth-order D and G operators are therefore (again only 
subscripts different from k, £, or m are explicitly shown). 


6u 6v 5w 


12A J u (k-2) " 8u (k-l) + 8u (k+l) 


- u 


(k+2) 


+ 12A | v U-2) ’ SV (il-l) + ^W+l) 


- v 


(£+ 2 ) 


+ 12A W (m-2) - 8w (m-l) + 8w (rH-l) 


- w 


(nH-2) 


div u + 0(A h ) 


(3.8) 


G(P) 


L 6 ^ 6 , 6 \ 

y e x 6x e y 6y e z 6z J 


e x 12A j P (k-2) ' 8P (k-l) + 8P (k+l) “ P (k+2)| 


+ e y 12A j P (£-2) ‘ 8P (£-1) + 8P (£+1) 


- P 


(£+ 2)1 


A 1 

+ e ,v A ■ P. -v - 8P, . . + 8P, , - . 

z 12A (m-2) (m-1) (tn+1) 


- P 


(m+2) | 


■ grad P + 0(A 4 ) 


(3.9) 
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where e , e , e are unit vectors in the x-, y-, and z-directions . 
x y z 


3.3 DG Operator 


4 ( 4 P ) 


DG(P) = T7- (3.10) 

Expanding one term of (3.10), using the fourth-order difference scheme, 


- 7 = 


(12A) 
+ 16P 


2 j P (k-4) " 16P (k-3) + 64P (k-2) 


(k-1) " 130P (k) + 16P (k+l) + 64P (k+2) 


" 16P (k+3) + P (k+4) | 


(3.11) 


If P is expanded in a discrete Fourier series, similar to (3.2), the 
Fourier Transform of (3.11) is identified as, 


/V 




, 14Ak l - 16 e 13Ak l + 64e 12Ak l + 16e 1Ak l 


- 130 + 16e- 1Ak l + 64e- i2Ak l - + e' 14 ^ P 


72A 


1 I 3 A. 

P 


2 j -65 + 16 cos(Ak^) + 64 cos(2Ak^) 


- 16 cosOAk^ + cos^A^) > P 


(3.12) 


(3.12) may be obtained directly from (3.4). Therefore, in Fourier space DG 
operator becomes 

A 

DG - -k| kj (3.13) 


Compare (3.12) to the following fourth order central differencing scheme, 

2 2 

which is a commonly used approximation to 3 /3x operator; 
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12A 


~ P (k-2) + 16P (k-l) " 30P (k) + 16P (k+l) 




15 - 16 cos(Ak^) + cos(2Ak^)> 


P 



P (k+2) | 

(3.14) 


(3.15) 


~ 2 2 ~2 
where k^ is defined by (3.15). k^, k| , and k^ are compared in Fig. 

3.2 for N = 16 . 


3.4 Transport Difference Operator 

Differencing the transport terms in the form of (2.28) will auto- 
matically conserve momentum in an inviscid flow. But the 
computation becomes unstable and the kinetic energy increases. This 
happens in real flows in spite of the dissipative nature of and 

the Leonard term. This non-linear instability, first reported by Phillips 
(1959) , arises because the momentum conservative form does not necessarily 
guarantee energy conservation, and truncation errors in the energy equa- 
tion are not negligible. 

Arakawa (1966) devised a differencing scheme that conserves both 
mean square vorticity and energy in two-dimensional calculations that use 
the vorticity and stream function as dependent variables. This is not 
useful in a three-dimensional flow. 

A fourth-order transport differencing scheme that does conserve 
energy and momentum was developed for the present work: 
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N 



j5_ 

6x 


(uf) 


' JK j ( “ £) <k+i) - <“ F) (k-i) 

+ “ ( *(k+D ' ? (k-i)> + F( “(k+i) - “(k-i) 5 } 

" 244 j (u£) (k+2) * <u£) (k-2) + u(£ (k+2) " £ (k-2)' 

+ £( “(k+2) ' "(k-2) 5 | (3 ' 16) • 


Again we only show subscripts that differs from k, Z, or m . For details 
see Appendix I. In the present work this was used for the terms 
3 (u^Uj)/3xj ; for the Leonard term a second-order version of (3.16) was 
used, 



_ 1 _ 

4A 



- (uf) 


(k-1) 


+ u(f 


(k+1) 


' f (k-l) )+ f(u (k+l) " u (k-l) } (3,17) 


The familiar second-order central difference approximation was used for 
V 3 in the Leonard term. 




(f 


(k+1) 



+ f (k-l) ) 


(3.18) 


For the residual stress terms, (3.17) was used. The strain-rates and 
vorticity were computed from the second-order central difference (3.6). 


3.5 Time Differencing 

A second order Adams-Bashf orth method was used for the time inte- 
gration. As shown by Lilly (1965), this method is very weakly unstable, 
but the total spurious computational production of kinetic energy is 
small . 

The Adams-Bashforth formula for u^ at time step n+1 is 
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-(n+1) 


n^ n) + At H^ n) - J + 0(At 3 ) (3.19) 


where is defined by (2.29). Note that this is an explicit scheme. 


3.6 Pressure Field Solution 

To study this problem in detail, let us rewrite (2.29) as 

S, k JL 

at ‘ h i " " 3x ± 

Again, continuity is 

3u 

JT = ° 


(3.20a) 


(3.20b) 


Taking the divergence of (3.20a) 


„2_ _J_ k _L 

v p = ax h i “ at a x 

i i 


(3.21) 


The usual computational procedure involves choosing the pressure 
field at the current time step such that continuity is satisfied at the 
next time step, i.e. so that the new flow field will be divergence free. 
This must be done very carefully. Let's look at three possibilities. 
These take advantage of Fourier transformation, for it is known that 
fast Fourier transforms provide an excellent way to solve the Poisson 
equation, at least in a rectangular domain. 

(a) Method 1 

The Fourier transform of (3.21) is 


2~ 
-k P 


(3.22) 


where g^ is the difference approximation to g^ and k 

are wave numbers in x-, y-, and z-directions . By inverse 


2 2 2 
k +k +k 
x y z 


k , k , and k 
x y z 
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transformation, P can be obtained. However, as will be shown shortly, 
this method does not give a next velocity field that is divergence- 
free in grid space. Hence this approach is unsatisfactory. 

(b) Method 2 

A second approach is to use the difference form of (3.21), 



(3.23) 


Then P can be obtained by Fourier transforming (3.23), 


*i*i * 


(3.24) 


As will be shown shortly, the pressure from this does not give a divergence- 
free field in grid space at the next time step. Hence, this method, which 
was used by Jain (1967), is also unsatisfactory. 

(c) Method 3 

The finite difference forms of (3.20a) and (3.20b) are 


where h^ is difference 
operator to (3.25a), 



Du^ = 0 

approximation to h^ . 


(3.25a) 

(3.25b) 

If we apply the D 


DG P 


if < B =i> 


+ Dh. 


(3.26) 


Then, taking the Fourier transform of (3.26), and using (3.13), 


-kjk^ P 


e' 

g l 


(3.27) 


Now let’s compare the three methods. The Fourier transform of 
(3.25a) is 
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(3.28) 


fiu ^ 

-sT'K 


-i k^ P 


To satisfy continuity in grid space, we want to have 

- 0 


5 

fit 


fi x i 


for a flow that has Du^ = 0 

0 at the first and next time steps 


^ w to start. From (3.5), this is equivalent 


Using (3.21), and sub- 


t0 

stituting P in (3.28) by P from (3.22), (3.24) or (3.27) and operating 
with D on the resulting equations, one obtains, 

for method 1: 


£ 




i 0 


(3.29) 


for method 2: 


£ (k iV 




t 0 


(3.30) 


for method 3: 


£ (k iV 


( s i - 0 


(3.31) 


The error introduced by method 1 and method 2 can be seen clearly by 

2 2 2 ~2 

observing the magnitude of the ratio k' /k or k' /k as illustrated 
in Fig. 3.2. Method 3 satisfies continuity at the next time step in 
grid space, and hence is chosen for pressure field solution here. 


3.7 Summary of Difference Equations 
Momentum equation: 


fiu^ 

fit 


fi 

fi Xj 

fiP 

fix. 


+ a a 

U,u,+ TT 


(u,u 4 ) - 


n 24 fix k fix k 


2v T S ij 


r __§l 

h i fix. 


(3.32) 
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I V 3 


Poisson equation for P: 


where 


P 

ij 


0 ) 


i 


DG(P) = D(h ± ) 


(3.33) 



(c $ A a ) 2 (2S ljS±j ) 1/2 

. . ,2 - . 1/2 

(C vV (u) i W i ) 



Smagorinsky model 
Vorticity model 


c„,c «= constant 
S v 
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CHAPTER IV 


DECAY OF ISOTROPIC TURBULENCE 


4.1 Problem Description 

Perhaps the most basic problem in turbulence is the decay of incom- 
pressible homogeneous isotropic turbulence. It is with this primitive 
turbulent flow that our study began. This flow was used to determine the 

value of the residual stress model constant (C or C ) , for use in sub- 

s v 

sequent calculations of other flows. It also provided a basic testing 
ground for the computational methods being developed in this program. 

The experimental grid turbulence data of Comte-Bellot and Corrsin 
(1971) were used as the "target" for these predictions. Such experiments 
closely approximate homogeneous isotropic turbulence, when viewed in a 
coordinate frame translating at the mean flow velocity. 


4.2 The Benchmark Experiment 

The pertinent information from Comte-Bellot and Corrsin' s (1971) 

experiments will be reviewed now. The wind tunnel test section, which 

has a slight secondary contraction to isotropize the turbulence, is 

sketched in Fig. 4.1 . The turbulence was generated by a biplane square 

rod grid with mesh size, M , of 5.08 cm. The free-stream air speed, U q , 

was 10 m/sec, giving grid mesh Reynolds number U M/v of- 34,000 . The 
2 2 2 ° 

streetwise (<u >) and transverse (<v > , <w >) turbulent energy components 
remained nearly equal to each other during the decay along the test sec- 
tion. These were closely fit by 



(4.1a) 


U 


<v 2 > 



(4.1b) 


Correlations, energy spectra and other quantities were measured using 
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hot wire anemometry at U Q t/M = 42, 98 and 171. The Reynolds number 

based on the Taylor microscale, = , was 71.6, 65.3 and 

60.7 at these points. 

4. 3 Mesh Size Selection 

The choice of the computational mesh size requires consideration of 
the turbulence spectrum (Fig. 4.3). The smallest scales that can be re- 
solved have wave number tt/A , where A is the mesh width. If there 
are N points in one direction, the largest scales that will be repre- 
sented in the computation have wave number 2 tt/(NA). N and A must be 
chosen such that the computation captures as much of the turbulence 
energy as possible. It is also desirable that the computation extend 
to the so-called inertial subrange (Tennekes and Lumley 1972). 

The mesh systems used were as follows: 

16^ mesh 

-3 

A = 1.5 cm, N = 16 , At = 6.25 x 10 sec 

3 

32 mesh 

-3 

A = 1.0 cm, N = 32 , At = 6.25 x 10 sec 

3 3 

The model constants were first evaluated using the 16 mesh; the 32 
calculation then verified that the constants are independent of mesh size. 
The corresponding Courant numbers were: 

2 —2 —2 —2 
q = <u > + <v > + <w > 


Nc 




q 3 /3 


At/Ax 


3 

16 mesh 

Nc <_ 0.06 

3 

32 mesh 

Nc < 0.1 
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A . A Initial Conditions 


We want to prescribe an initial profile that has the proper energy 
content and spectrum, and is isotropic. A technique for generating a 
random field that meets these conditions, and also satisfies continuity, 
was developed for this purpose. Since the computation will treat the 
filtered field, we matched the initial filtered spectrum and not the 
Initial measured spectrum. 

We generated initial filtered field, u. (x) , by first establishing 
its discrete Fourier transform, u. (k ) , 



u(x) = ELE (4.2) 

N 

n = - 2 


Here k is the wave number vector defined by (k ) , = rnr n , where n. 

— n ' n i i i i 

2 

is an integer ranging from -1/2N to 1/2N-1 for an N mesh system. Note 

that the maximum wave number is k = ir/ A . If u is discretized at 

max — 

N points, then the Fourier transform u can only be evaluated at N 
discrete wave numbers, and that is why the summation must have non-sym- 
metric limits. 

The commonly used fast Fourier transform requires N to be 2 m , 

where m is an integer (see Cochran et al. 1967). Physically N has 

to be large enough so that wave-number spectra can be treated as smooth 

3 3 

functions. As will be shown later, 16 or 32 mesh systems gave fairly 
smooth three-dimensional energy spectra. 

Now, we can approximate the spectrum function, of the filtered 

field (see Tennekes and Lumley 1972) as, 

$,.,00 ~ <u, (k )u *(k )> (A. 3) 

ij n l - n j n 

where < > denotes an average over an esemble of experiments or, alter- 
natively, over a spherical shell in k-space with radius k^ (see Section 
A. 6). 
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The filtered 3-D energy spectrum, E , is given by 

E(k) = 27Tk 2 $ il (k) (A. 4) 

E(k)dk is the energy content of a differentially thick spherical shell 
in wave-number space. Using (4.3), E(k) is approximated by 

E(k ) « 2irk 2 < u. (k )u *(k )> (4.5) 

n n i n i n 

To establish the initial field we need to fix the Fourier amplitudes 

A A A 

u. (k ) . Equation (4.5) was used to fix u. (k )u. *(k ) for each k . 
in l n i n 

The vector components u^ were chosen by a technique described below. 

/\ 

To get u. to satisfy continuity in grid space, the real and imag- 

inary parts of the transformed velocity vector, ij , must be perpendicular 

to the modified wave-number vector, k* (see Equation (3.5)). In the 

3 

actual computation, we have N points in k-space, and, for any k , 

A 

k' can be obtained by (3.4). Then, u has to be selected on a perpen- 
dicular plane to k' in k-space. 

/s 

To ensure statistical isotropy, the real and imaginary parts of 31 
must be chosen randomly. First we picked a unit vector A , perpendicular 
to k’ , by turning a random angle, 4> , from a reference frame (Fig. 

4.2). Here <j> was selected with uniform probability over the interval 
0 to 2ir . We then repeated to get a second random unit vector JS , 

A 

also perpendicular to k* . The real part of u was made proportional 
to the vector A and the imaginary to IS , and hence continuity was sat- 
isfied. We still needed to fix the relative magnitudes of the real and 

A 

imaginary parts of u(k) , which we did by a random choice of an angle, 

0 . Then we defined a and b by 

a = cos9 , b = sin 6 (4.6) 


Finally, we set 


■ Msh )0 * < aA j + lbB J > 


(4.7) 
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Now, by inverse transforming , we obtained , 
real and will be real if the Fourier transform satisfies 


which must be 


u(-k ) = u (k ) 


(4.8) 


In essence, the imaginary contribution for each negative k 

XI 

cancels the imaginary contribution of the same positive k 


n 


exactly 
Hence, we 


only needed to generate u. by (4.7) for the upper half of the k-space, 


i.e. for 


N 

0 < n < | - 1 


However, this won't fix u for 


N 

n i 2 


Moreover, if u is not zero, the velocity field will have an imaginary 
part (see (4.2)). If instead we wrote (4.2) as 


u ± (x) 


N/2 

E E E u i ( ^n )e 

-N/2 


ik 


then u. would be real. However, then we could not take advantage of 

. As a practical solution to this dilemma 


the FFT routine to invert 
we set 

-1/2N . Then (4.2) is essentially the same as 


Uj equal to zero for the wave numbers corresponding to n^ = 


(f - 1) 

u(x) = EEE u(k )e 

-(§ - 1 ) 


ik 


and u^ will be real, and an FFT routine may be used. The resulting 
energy spectrum was therefore slightly low at the highest wave number. 
However, the effect of this discrepancy was insignificant and became 
invisible after a few time steps in the computation. 

We remark that the field generated by this procedure is quite iso- 
tropic. However, as will be shown it has zero skewness, whereas real 
turbulence has a non-zero skewness. As will be seen, this condition 
corrected itself in only a few time steps. 
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4. 5 Boundary Conditions 

The computational problem can only extend over a part of the experi- 
mental region. To get around this difficulty we have used "periodic" 
boundary conditions, which are of course not really correct. However, 
if the computational grid system extends over a distance large compared 
to the scale of the energy- containing motions, the periodic boundary con- 
ditions should not introduce appreciable error. The periodic boundary 
conditions have a great advantage in dealing with the Poisson equation 
for the pressure by fast Fourier transform (FFT) methods. 

4.6 Extraction of Statistical Information from the Computation 

The statistical quantities of interest are averages over ensembles 
of experiments. Since we made only one computational realization in 
each case, the statistical quantities had to be inferred from appropriate 
ergodic hypotheses. 

In physical space the ensemble average < > was replaced by an 
average over the flow field. This was done by taking a mean value over 

3 

N mesh points, i.e. 


N 

<f (x) > - Jjj EEL *(*.».-> w-9) 

k,£,m=l 

The differencing schemes described in Chapter III were used to calculate 
these quantities. 

In wave number space, the ensemble average was replaced by an aver- 

3 

age over a shell in k-space ("shell average"). Since we have only N 
discrete points in k-space, the < > average was made by taking a mean 
value over the points between the two shells with radius (k-l/2Ak) and 
(k+l/2Ak). 

A 

To get the filtered spectrum, E(k) , the transformed velocity, u , 

A — * 

was obtained by FFT (see 4.2). Then, <u^(k^)u^ Qi n ) > was calculated 
by shell-averaging. The choice of the band width, Ak , is somewhat 
arbitrary and was set to be 0.1 cm ^ here. Then, from (4.5), 
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k+Ak 


— 2ttk X"' — — * 

E < k > 2^ u i ( V u i 


— n 


(A. 10) 


k =k-Ak 
n 


-1 


where is the number of points between the two shells with radius 

(k-l/2Ak) and (k+l/2Ak) . The resulting spectra, evaluated at 0.1 cm 
wave-number intervals, were smooth enough to be represented by continuous 
curves as shown by solid lines in Fig. A. A, A. 5, and A. 7 . The filtered 
spectrum, E(k) , was then compared to the filtered experimental spectrum, 
which we obtained using (2.23). 


A. 7 Selection of the Averaging Scale 

Considerable thought was given to the choice of the averaging scale 
A^ . Our failures are as important as our successes, and both will now 
be discussed. 

Consider first the computation with A^ = 0 . This zero averaging 
length is equivalent to the unfiltered calculations used in laminar flow, 
and implies that we are trying to resolve the complete spectrum by a fin- 
ite difference method. The Leonard term in this case is equal to zero, 

1 . e . 


U i U J 


U i U j 


The unfiltered initial energy spectrum is plotted in Fig. A. 3 . The 

3 3 

amounts of unfiltered energy for the 16 mesh and 32 mesh systems are 

also shown. Figure A. A shows the computation for a value of C that 

s 

gives the proper rate of energy decay. Note that for k > 1/ 2k ^^^ the 
spectrum is distorted considerably at tU^/M = 86.5 and become worse as 
time increases. 

In an instantaneously fluctuating field, higher derivatives are not 
small and the convergence of the Taylor series is expected to be slow. 

Use of A^ = 0 and the consequent exclusion of the Leonard term caused 
much distortion of the spectrum, i.e. aliasing error. Indeed, the fin- 
ite difference method with N mesh points in one direction can only re- 
solve the unfiltered field up to k = it/ (2A) = k /2 , which is a half 

niflx 

the maximum wave number in one direction (see Orszag 1969, Orszag and 
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Israelii 1974). Judging from Fig. 4.4, in a computation without filter- 
ing the non-linear interactions transfer too much energy from large to 
small scales. This excess up-scale energy transfer could be somewhat re- 
duced by using a larger coefficient in the residual stress model. How- 
ever, this brings an unreasonably high energy decay rate. We conclude 

that one should never use A, = 0 in a turbulence simulation; filter- 

A 

ing is essential. 

Let's look next at what happens when the averaging scale A^ is 

equal to the mesh scale A . Figure 4.5 shows a computation with a value 

of C that gives the proper energy decay. Note that the errors in the 
s 

predictions of the filtered spectrum are significant at high wave numbers. 

Consider now the computations with A^ = 2A , shown in Fig. 4.6, run 
for values of C g and that give the proper energy decay. The pre- 

dictions of the filtered spectrum for both residual stress models are 

3 

remarkably accurate, even in the coarse 16 calculation! 

To Investigate the effect of the Leonard term separately from the 
filtering, an additional calculation with A^ = 2A was run with the 
vorticity model, excluding the Leonard terms (Fig. 4.7). The prediction 
is poor on high wave number side. Evidently the Leonard terms assists 
in removal of energy from high wave numbers. We conclude that good re- 
sults will be obtained with A A = 2A , and that the Leonard terms must 
be included. 

4.8 Selection of C and C 
s v_ 

An analytical way of determining the residual stress model constants, 
C g or , is not known. Lilly (1966) estimated C g = 0.2 using sev- 
eral ad-hoc assumptions. Later workers (Deardorff 1971, Fox and Deardorff 
1972) calibrated this constant to get the best computational results. In 

these cases, the required C was between 0.10 and 0.22 . 

8 3 

In the present study a series of 16 mesh calculations were run with 
different values of each constant, and values selected that gave the best 
prediction for the filtered rate of energy decay as judged by consideration 
of the slope of the curve (Fig. 4.10). The constants obtained were as 
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C = 0.254 

v 

Figure 4.9 shows the sensitivity of the predicted filtered energy decay 
to . Figure 4.10 shows the excellent agreement of the energy history 

with the data for the final constants. 

3 

Figure 4.11 shows the energy decay rate from the 32 calculation 
with these same constants. The spectral results are shown in Fig. 4.8 . 
The excellent agreement with data confirms that the model constants do 
not vary with the mesh size, at least in the range covered. 

In comparing Fig. 4.10 and 4.11, it must be remembered that these 
are the filtered energies, which are different in these two mesh systems. 
Because of the discrete Fourier approximation, not all the turbulent 
energy is captured (see Fig. 4.3). Filtering improves the situation, be- 
cause less energy is omitted from the filtered field at high wave number. 
However, the energy in the discrete approximation to the filtered field 
was still less than that in the filtered experimental field. To facili- 
tate selection of the constants, the filtered experimental history was 
shifted as shown in Fig. 4.10 and 4.11 . 

4. 9 Energetics of the Filtered Field 

Multiplying (2.28) by u^ , and taking an ensemble average, and as- 
suming homogeneity, one finds 

- - £ ■- te H + e i> 

2 — — 

where q = <u^ u^> . 

The dissipation e is seen to have two parts, a part representing 
transfer to the residual field. 


follows ; 


0.206 


The three digits are not meant to imply accuracy. We actually ran 
with 

(2C g ) 2 =0.17 , (2C y ) 2 =0.26 . 
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(4.11b) 


e R - - <u i 3^ < 2V T V" 


and a Leonard term part. 


c , - 3 , A „2 - - . 

£ l ■ + <u i v u i V 

j J 


(4.11c) 


To see the relative contributions of the Leonard term and the re- 
sidual scale motions to the energy decay rate, e,/e and £_,/£ are 

L o K o 

shown in Fig. 4.14. Here e is the sum of £_ and £„ at (U t/M - 

O L K O 

3.5) = 42 where the computation started. £ is much smaller than 

L 

Leonard estimated (1973) . As shown by Leonard (1975) , £ takes energy 

mostly from the large wave number side thus preventing the damming up of 
energy in the smaller eddies. 


4.10 Other Aspects 

No significant difference is observed between Smagorinsky and vor- 
ticity models. However, some differences are expected in future applica- 
tions to unbounded flow problems with turbulent and non-turbulent regions. 

The skewness, which is a measure of vorticity production in the 
energy cascade process, is shown in Fig. 4.12 and 4.13. Since the initial 
field is randomly generated, the skewness is zero initially, but quickly 

3 

adjusts to essentially a constant value. For the 16 calculation the 
value is clearly too low (the experimental skewness is about-0.4). For 

3 

the 32 calculation the skewness seems slightly high. 

We have emphasized the need for a fourth order differencing scheme, 
and wonder why others have been able to do so well with second order 
schemes. The reason may be that the second order difference form of the 
advection term implicitly includes Leonard-like second order truncation 
terms and thus the Leonard term is partially taken care of by the trunca- 
tion. If a fourth order scheme is to be used, the Leonard terms should 
be included explicitly. We have seen that they are important, particularly 
at the high wave numbers. We conclude that, for a grid calculation of 
the type run here, the best results will be given by the fourth-order 
difference scheme that incorporates the Leonard terms. 
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A question arose as to the behavior of the vorticity under the dif- 
ference scheme used here. A two-dimensional irrotational flow was input, 
was set equal to zero, and two time steps were taken. The vorticity 
remained exactly zero, indicating that, at least in a two-dimensional 
flow, the differencing scheme will not produce unwanted vorticity. This 
aspect of the computation should receive further study in the future. 

4.11 Computational Details 

The calculations described above were executed on the CDC 7600 at 
Lawrence Berkeley Laboratory, using programs written in FORTRAN (Appendix 
II). The total storage requirements (octal) for 60 bit words were as 
follows : 

3 

16 calculation 

Large Core Memory: 230,360 

Field Length (Small Core) required to load: 121,200 

3 

32 calculation 

Large Core Memory: 1,100,234 

Field Length (Small Core) required to load: 121,200 

The computer time per computational step was approximately as follows 
3 

16 calculation : CPU time = 3 sec 

3 

32 calculation : CPU time = 20 sec 

The calculation program was carefully checked before these production 
runs. To check each term in the difference equation (3.32) and (3.33), 
we imposed systematically artificial flow fields. For the terms involving 

— g 

first derivatives of velocities such as S.. , a). , v_ , -= — (u. u. ) and 

i j i T ox^ i j 

, the following linear velocity field was used: 

u = x + 2y + 3z 

v = 4x + 5y + 6z 

w = 7x + 8y + 9z 
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Then the computed results were compared to the exact values. For the 
terms with second derivatives, the following quadratic expressions were 
used: 


u = 

2 

X 

+ 2y 2 

+ 

3z 


, 2 

2 



V = 

4x 

+ 5y Z 

+ 

6z 

_ 

2 

2 



w = 

7x 

+ 8y 

+ 

9z 


Then, at randomly picked mesh points, the computer results for the ad- 
vection terms, the Leonard terms, fi^ , and D(fi^) were compared to the 
exact values obtained analytically. 

For the Poisson solver, a sinusoidal pressure field was used to 
generate DG(p) , then the computer results were checked against the im- 
posed pressure field. The initial field was generated as described in 
Section 4.4 and two time steps were advanced. The subsequent results 
provided a testing ground for time stepping, the maintenance of a di- 
vergence-free velocity field, the overall sequence of computing, and 
input, output, tape handling, and data reduction routines. 

The computer program is given in Appendix II. 
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CHAPTER V 


DISTORTION OF HOMOGENEOUS TURBULENCE BY 
IRROTATIONAL PLANE STRAIN 


5.1 Problem Description 

Shear may be viewed as a combination of pure strain and rotation. 
Therefore, a basic problem is that of homogeneous turbulence acted upon 
by an imposed uniform homogeneous irrotational strain. Tucker and Reynolds 
(1968) approximated such a flow experimentally by passing grid-generated 
turbulence through a passage designed to produce uniform strain in a 
coordinate system translating with the mean flow velocity. This experi- 
mental flow approximates the problem of box turbulence with a constant 
rate of strain, shown in Fig. 5.1b. 

In this chapter we discuss the computation of an idealized homo- 
geneous flow with irrotational pure strain, comparable to the Tucker- 
Reynolds laterally strained flow. In addition, we treat the return to 
isotropy following the removal of strain, which roughly corresponds to 
the experiment in the uniform channel downstream of the straining section. 

Tucker and Reynolds did not measure the energy spectrum and hence we 
cannot make an exact comparison with their data. However, for a qualita- 
tive comparison, the initial turbulent intensities in the computation 
were set to be equal to the experimental values at the beginning of the 
strained section. Two cases were run. The first case was run with 
approximately the same initial anisotropy as the experiments. However, 
there are problems in that the anisotropic field so generated had improper 
shearing stresses. Therefore, a second calculation was made with an 
initially isotropic field, and this flow has been used to study the effects 
of pure strain on homogeneous turbulence. 

The initial field for the computation was based on an energy spectrum 
similar to that used in Chapter IV. However, the Tucker-Reynolds initial 
energy level was much higher than the energy in the grid flow studied in 
Chapter IV. To adjust the energy, the amplitude of the Fourier coefficients 
were multiplied by a constant. The initial one-dimensional energy spectra 
are shown in Figs. 5.11 and. 5. 12 by solid curves. 
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The strain-rate used In the calculation was different from that of 
the Tucker-Reynolds experiments. In an initial calculation their strain 
rate was used, but the energy-decay rate in the relaxation section did 
not match their experiments. The difference was attributed to differences 
between their (unmeasured) spectrum and that used to start the calculation. 
We first computed the energy-decay rate in the absence of strain as shown 
by the solid lines in Fig. 5.5 and 5.8. Then, the strain -rate was deter- 
mined to get the same total strain as in the Tucker-Reynolds experiments 
at a point in the flow that would have the same energy in the absence of 
strain. The final calculations were performed using this strain rate. 
Therefore, the calculation should be regarded as a "Tucker-Reyno lds-like" 
flow, and not as a simulation of their flow. 


5.2 Governing Equations 

To handle the imposed mean strain, we express the local velocity 
and pressure field as* 


u i (x,t) » U^x) + u^'(x,t) (5.1a) 

p(x,t) - P(x) + p"(x,t) (5.1b) 


where 

u ± “ (U,v,0) - (rx,-r y ,0) (5.2a) 

P - - Y pr 2 (x 2 + y 2 ) (5.2b) 

T is the constant strain-rate. With this decomposition, the Navier- 
Stokes equations for incompressible flow become 


_ 3 _ 

dt 




(5.3a) 


*Note that we place the strain in the pl ane > while Tucker and 

Reynolds placed it in the x^-x^ plane. 
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(U ± +u J) 


0 


(5.3b) 


3 

3x, 


2 

Now, using the definition (5.2), and noting that 3 U^/3x^3x^ *» 0 , 

this reduces to 




0 


(5.4a) 

(5.4b) 


These are the equations that will solve by the methods presented pre- 
viously. The only modification (compare (2.5)) comes from the third 
term on the left; this term may be regarded as a forcing function due 
to the mean strain. 

Now, we express each variable quantity, f" , as 


where 


f" - f + f' 

f(x) = /g< 2 -x') f"(x') dx* 


(5.5a) 

(5.5b) 


Note that f is now the filtered f" field. The U^u^' term in (5.4a) 
is filtered using the method described in chapter II, giving 


u i u }'<V 


r 

r 


3U i (X n> 

1 

■ y g(x,2^) 

b 


+ ^ ° 
3X k 

( V*lco ) l 

o ■* 

[W + 

3u . 

<*„> 

(x k’ x ko ) + 

1 3 “t (x ' 

2 3x k 3 £ 

+ . . .u! dx 





■ U i u j + 24 

V 2 

Vj 

+ 0(A*) 



(x k" x ko ) (x r x iio ) 


(5.6) 
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-The model for the residual stress, , is taken as 


l ij ° 3 R U 6 ij “ V T 2S ij 


(5.8) 


where now is the total strain-rate. 


3 ij - i|^- (u i ^i )+ 4 ( vyj 


(5.9a) 


Since we expect the strain-rate to be dominated by small scales, the 
eddy viscosity, v T , was evaluated using the vorticity model, 


V T ‘ (C vV 2 <“l“i )1/2 


(5.9b) 


o) ■ curl u 


(5.9c) 


Since the imposed flow is irrotational, v T is based on the total RMS 
vorticity. 

Then, filtering (5. A), and again neglecting the viscous term, the 
following equations are obtained (compare (2.28)). 


_a_ * 4 

at “ " axj ( u i u j + 2 A 


v 2 (u ± Uj ) 


2 v t (s ij + ^> 


9P 

Sx, 


+ F. 


(5.10a) 



0 


where 


and 



(5.10b) 

(5.11a) 

(5.11b) 
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F^ are the terms through which the major effects of the mean strain come 
in. These are 


F 1 • ^ [vj + Vi + 57 ’ 2 <WYi>] 

Note the appearance of a Leonard correction term. 

For the computation, the difference form of (5.10) is used, 

k2 o \ 

M 


(5.11c) 


6u i 6 i- - A A 

"""fit" ° -6^{ u i u j + 24 


2v_(S..+, 


v\ ij T 1J 


- G(P) + F. 


(5.12) 


where F^ is the difference form of F^ and 

s . u 6 A + s h 

S ij 2 ^ 6xj &x ± 


(5.13) 


The exact expressions for ^ ^ were used. As before, the Poisson 
equation for P is obtained by operating with D on (5.10a), 


DG(P) 


D(h ± ) + D^) - D(u ± ) 


(5.14) 


where 


_fi 

fix 


j 


u i u j + 


_A 

24 fix, fi 


5 -- 


V*k 


u i u j 


- 2v x (S ij +J lj ) J (5.15) 


The space differencing and the time advancing schemes are the same as 
those explained in Chapter III. Periodic boundary conditions in all 
three directions were imposed, and the same solution procedure as 
described in Chapter IV was applied. 


5.3 Anisotropic Initial Condition 

Anisotropy in grid generated turbulence is not negligible in many 
experiments (e.g. Grant and Nisbet 0-957) , Uberoi (1963), Tucker and 
Reynolds (1968)). Therefore, to make the initial condition reasonably 


39 



close to the experiments, It was felt desirable to generate an Initial 
field with anisotropy. A method for this will now be described. 

Suppose that the u and v components of turbulence energy are 
equal while the w component is different, 


<U 2 > “ 

<V 2 > 

(5.16a) 

—2 

<w > = 

(1+R) <u 2 > 

(5.16b) 


where < > denotes an average over an ensemble of experiments. For 
the Tucker-Reynolds experiments , w is in the mean flow direction and 
R - 0.45. Now let us decompose w into its isotropic part, w^. , and 
anisotropic part, w^ , such that 

— 2 — 2—2 

< u > * < v >-<Wj> (5.17a) 

and 

w * w T + w. (5.17b) 

1 A 

If we assume that the isotropic part can be generated by the method in 
Chapter IV, then, for continuity to be satisfied, 

D (w A ) = 0 (5.18) 

This is a crude assumption. However, unless we know more about the 
initial turbulence structure, this is perhaps the best we can do. Now 
(5.18) in Fourier transformed space can be written as, 

A 

w A k^ = 0 (5.19) 

where ~ denotes a Fourier transform and kl is defined by (3.4). 

A J 

Therefore, w can have a non-zero values only when k^ = 0 . Then, 
following the same procedure discussed in Section 4.4, we get 
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Ca+ib) 


(5.20) 


W A (k l ,k 2 ,0) 



where 


~2 


q 


A/N A A A A A A 

uu* + w* + W T W* + w.w* 
II A A 


and a ■ cos 0 , and b ® sin 0 are obtained from a random angle, 0 , 
with uniform probability from 0 to 2ir . 

The initial condition for the first run was generated by this pro- 
cedure and, for an input R of 0.45 , the generated field emerged with 
R = 0.43 . This field had shearing stresses not present in the actual 
flow, and so the second run was made using an isotropic initial field 
generated by the method described in Chapter IV. Further studies in 
Section 5.5 and 5.6 are based on the second run. Both runs are reported 
for completeness. 

I 

5.4 Results 

The results of the following two cases are presented. 


(1) Anisotropic 


initial field 

-2 -2 
<w > <w > 

-2 ~ -2 
<u > <v > 


1.43 


(2) Isotropic initial field 

-2 -2 -2 

<u > = <v > = <w > 


For both cases, the mean stream speed was taken as w 


T •» 1.457/sec., A = 0.59 in and At 


5.36x10 3 sec. 


240 in/sec, 
This corresponds 


to a Courant number 

V q2 / 3 + U nax At/,A - 0,15 ' The 16^ mesh system 
was run using the vorticity model, with = 0.206 as obtained from 
the isotropic decay studies. 
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For convenience of reader, three of the Tucker-Reynolds measurements 

for the turbulent intensities, <u? > /w* the turbulent energy ratios, 

2 2 ^ ® . _2 o 2 2 

<u^ > /q , and the structural parameter, = {<v > - <u >}/{<v > + <u >} , 

are replotted in Figs. 5.2, 5.3 and 5.4 For the case 1 , the same quanti- 
ties are plotted in Figs. 5.5, 5.6 and 5.7. Here, for comparison, the 
time in computed results, t , are converted into the downstream distance, 

Z , by Z « W^t . The behavior of these are comparable to Tucker and 
Reynolds results. However, when the strain is turned off, the rate of 
return to isotropy is much slower than in the experiments. It is inter- 
esting that this is consistent with the feelings of some workers that the 
return to isotropy in the Tucker -Reynolds flow is too rapid (Reynolds 
1975) , perhaps because of defects in the experimental simulation of 
homogeneity. The same quantities for case 2 are shown in Figs. 5.8, 

5.9 and 5.10. 

One-dimensional energy spectra were obtained in a similar manner 
to that described in Chapter IV (see Tennekes and Lumley (1972)). The 
only difference is that the shell average in Chapter IV is replaced by 
a plane average In wave-number space, i.e. 

-ip? "<v “*<v < 5 - 2i > 

N k 2 k 3 

Here, the notations are the same as before. These spectra were computed 
at three different times or downstream locations (Figs. 5.11 and 5.12). 

At zeroth time step, » ^ 22 ^ 2 ^ anc * ^33 ^ 3 ) are almost identical, 

as they should be. By the end of the straining period (75th time step), 

E^ and E 22 have become quite different. is flatter on the large 

eddy side while both small-scale spectra are nearly the same. Over the 
last period, in which there is no strain, the spectra approached one 
another very slowly, as seen by the spectra at the 125th time step. 

The calculations were run on a CDC 7600, and required approximately 
7 minutes for each case. Storage requirements were similar to those for 
Isotropic decay. 
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5.5 Energetics of the Filtered Field 

Multiplying (5.10) by u^ , taking the ensemble average, and 
assuming homogeneity, one finds (compare (4.11)) 


G* + Cl 


(5.22) 


where 


q 

<P 

(?L 


<u iV 


-<u.u > ~ 

i j 


3U 


(5.23a) 

(5.23b) 


‘j 


2* J <u i 3 ^ V2 <V j )> + <u i 9 ^ V2 (U j U i )> > (5 ‘ 23c) 


e R and are the dissipation terms discussed in Section 4.9. (p and 

are the production terms. Note the appearance of a Leonard production 
term, (p L » that comes from the non-linear interaction between the mean 
strain and the filtered field. 

The computed behavior of these terms is shown in Fig. 5.13. Note 
that <Pl contributes significantly, particularly where the anisotropy 
is large near the end of straining period. 


5.6 Assessment of Turbulence Closure Models 

Turbulence computation of the conventionally averaged (ensemble or 
space) quantities has been based on some ad hoc closure models with a 
number of adjustable constants. As pointed out by Reynolds (1974a, b), 
more systematic approaches are desirable for generalized turbulence 
models. Even though laboratory experiments provide actual quantities, 
experiments are limited because important properties like the pressure- 
strain correlations are difficult to measure directly. On the other 
hand, computer-generated experiments provide a vast amount of data on the 
flow field, and hence the numerical experiments can be used to study 
the closure models. We have attempted to study the pressure-strain terms 
and other statistical quantities using the present computation. Even 
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though the 16x16x16 mesh calculation gives good results for the energy 
components, as will be shown shortly we cannot use the computer generated 
field to compute the pressure-strain term directly, at least not in the 
present calculation. 

The exact Reynolds stress equations for homogeneous flow without 
mean deformation are 


dR 


where 


JA 

dt 


Ij 


ij 


T - D 

ij ij 

3u 


V ax j Sx i/ 


3u. 3u. 


2\> < -r-^- -r— L > 

3 *k 3 *k 


(5.24) 


Here u^ denotes the turbulent components of velocity and p is the 
turbulent pressure divided by the fluid density. 

The pressure-strain term T^ is responsible for the return to 
isotropy following removal of strain. The modeling of T^ has been 
the subject of much discussion. Since no direct measurement of this 
term is known, we tried to estimate this term using the present compu- 
tation in the retum-to-isotropy portion of the computation. 

was obtained from 


The computed pressure-strain term, (T^ ) 


( V= 


< P 


,3u i + ^ 

V 9x j 3x J 


(5.25) 


Here < > denotes an average over the flow field. The fourth order 
central differencing scheme was used for the 5/5x. terms. 


j 

For comparison T^ was obtained a second way. Using D 


ij 


2/3 e 


(see Reynolds (1974b), (1975)), T^ was obtained from the computed R 


history as 


ij 


( Vr “ £ <u iV + l £S u 


(5.26) 
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e was estimated from the energy equation, which in the absence of strain 
is 


d <u^u^/2> 
dt 


(5.27) 


e agreed well with £ + £ , computed directly. The time derivatives 

R Li 

were approximated by a second order central differencing formula. 

Finally, we predicted by Rotta's (1951) model for T^ in the 

absence of mean strain, 


( Vm 


-A £ b . . 
o ij 


(5.28) 


where A is a constant and 
o 


’ij 


<u iV iii 

2 " 3 

q 


We used A q “ 2.5 , as suggested by Reynolds (1975). 

These three results are shown in Fig. 5.14. It appears that (T .) 

ij c 

is quantitatively poor. This is attributed to the coarseness of the 

3 

16 mesh . We conclude that T^ undoubtedly contains some Leonard-like 

terms, and this must account for the difference between (T J .) and 

ij c 


<Vr • Whlle T 1J 


can be estimated from the R^ equation a la 


(5.26), it cannot be computed directly from the calculated field with 

3 

such a coarse grid. A repeat of this work using a 32 grid is recommended. 
This should be accompanied by a careful analysis of the Leonard terms 


arising in the R^ equations. 
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CHAPTER VI 


CONCLUSIONS AND RECOMMENDATIONS 

In this thesis we have developed the basic approach to computation 
of three-dimensional time-dependent turbulent flows. We have seen that, 
with a modest 16x16x16 mesh and a residual scale stress model, many 
interesting features of experimental flows can be computed. Work remains 
to be done in the development of better approaches for using this type 
of computation to assess turbulence model equations, and to extend the 
procedure to other flows. 

It would be informative to study the effect of the initial spectrum 
on the rate of return to isotropy. This might be done by removing the 
strain at different points along the straining portion of the Tucker- 
Reynolds flow, allowing isotropy to be restored from points with different 
spectra. 

In extending the method to other interesting flows, problems to be 
resolved include the handling of non-periodic boundary conditions, solid 
boundaries, and free boundaries connecting the region of computation to 
irrotational flow outside. 

One useful problem to study would be the case of homogeneous tur- 
bulence near a wall, without mean shear, for which some experimental data 
exist (Uzkan and Reynolds (1967)), and would be the diffusion of tur- 
bulence into a non-turbulent region, again without shear. It is recommended 
that experience with these simple problems be gained before a more com- 
plex flow is attempted. 

When one moves on to handle flows like jets, wakes, and mixing regions, 
it should be possible to take advantage of the fact that the flow outside 
of the superlayer is irrotational, and to use the exact solution for 
unsteady irrotational flow to extend the calculation to infinity out 
beyond the mesh. Care must be taken that the numerical scheme does 
not produce vorticity in an irrotational flow; the diffusion of vorticity 
by will also have to be handled in a way that prevents its diffusion 

into the irrotational external field. 
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Eventually It may be possible to treat practical flows, such as 
boundary layers, wakes, combustion, etc.; by these methods. But much 
more effort should first be devoted to fully understanding the nuances 
of the residual scale models, grid-schemes, differencing schemes, filters, 
etc. that are the bases for this type of numerical simulation. 
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Fig. 2.1. Gaussian Filter 

G i ■ £ exp{ -6(x i -x|) 2 M 2 } 
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Fig. 3.1. Comparison of Modified Wave Numbers 
(see Eqn. 3.4 and 3.7) 
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2 


Fig- 3*2* Comparison of Div(Grad) Operator 
= exact 

~~2 

k 2 “ fourth order (see Eqn. 3.15) 
2 

k' - fourth order (see Eqn. 3.12) 


GRID 
(mesh M) 




Uq = lOm/sec 


U = I.27U 0 


f— l8M-«-| 


Fig. 4.1. Sketch of Wind Tunnel Test Section for a 
Generation of Isotropic Flow. 
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Fig. 4.7. Filtered Energy Spectra-Effect of Leonard Term 
16x16x16 Mesh: Vorticity Model: 

A, = 2 A ; A ■ 1.5 cm 
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Fig. 4. 


8. Filtered Energy Spectra 

32x32x32 Mesh: Smagorinsky Model 

A. = 2A ; A = 1.0 cm 
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Fig. 4.10. Decay of Mean Square Filtered Velocity 

16x16x16 Mesh: = 2A ; A = 1.5 cm 

< > : Average Over all Space 
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Fig. 4.11. 
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Decay of Mean Square Filtered Velocity. 
32x32x32 Mesh: Smagorinsky Model: 


- 2A 
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Fig. 4.13 Comparison of Skewness, S . (see Fig. 4.12 for 
Definition of S ): Smagorinsky Model: 

A a - 2A 
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Fig. 4.14. Component 

16x16x16 
Vorticity 
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Fig. 5.1a. Schematic of Wind Tunnel Producing Constant 
Rate of Strain in x-y Plane 


V = -TY 



Fig. 5.1b. Equivalent Representation of the Plane Strain 
in a Box T ° 1.457 
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Fig. 5.2. Tucker-Reynolds Flow. 

Turbulent Intensities. The solid lines are for 
decaying turbulence in the absence of strain. 

Z = Downstream distance from the beginning of 
strained section. 
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Fig. 5.3. Tucker-Reynolds Flow 

Turbulent Energy Ratios 
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Fig. 5.4. Tucker-Reynolds Flow 

Change in Structural Parameter, 

-2 -2 

v <v > - <u > 

K 1 = -2 -2 

<v > + <u > 
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(Z + 10) (in) 

Fig. 5.5. Simulation of the Tucker-Reynolds Flow 
Turbulent Intensities. The solid lines 
are for decaying turbulence in the absence 
of strain. (Anisotropic Initial Condition) 
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Fig. 5. 
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6. Simulation of the Tucker-Reynolds Flow 
Turbulent Energy Ratios under the Plane 
Strain and the Return to Isotropy in 
Parallel Flow. (Anisotropic Initial 
Condition) 
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Fig. 5.7. Simulation of the Tucker-Reynolds Flow 
Change in Structural Parameter. 
(Anisotropic Initial Condition) 
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Fig. 5.8. Simulation of the Tucker-Reynolds Flow 
Turbulent Intensities. The solid lines 
are for decaying turbulence in the absence 
of strain. (Isotropic Initial Condition) 
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Fig. 5.9. Simulation of the Tucker-Reynolds Flow 

Turbulent Energy Ratios Under the Plane Strain 
and the Return to Isotropy in Parallel Flow. 
(Isotropic Initial Condition) 
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5.10. Simulation of the Tucker-Reynolds Flow 
Change in Structural Parameter, 
(Isotropic Initial Condition) x 







Fig. 5.11a. One-Dimensional Energy Spectra 
(Isotropic Initial Condition) 
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Fig. 5.11b. One-Dimensional Energy Spectra 
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Fig. 5.11c. One-Dimensional Energy Spectra 
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Fig. 5.12a. One-Dimensional Energy Spectra at 

Three Different Downstream Locations 
(Isotropic Initial Condition). 
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Fig. 5.12c. One-Dimensional Energy Spectra 











APPENDIX I 


ON THE FOURTH ORDER CONSERVATIVE SPACE DIFFERENCING SCHEME 


To explain the difference formula used here in detail, consider the 
following equations. 


3u 


9r + ^ ( vy ■ ° 

3u i 

" 0 


(A. 1) 
(A. 2) 


u ± ♦ (A.l) 


3u i . 3 

u i IT + U 1 3^ Vj 


3 / u i u lV 1 3 x Vl ^1 
3? \~T-J* 2 SJ Vl u J + 2 3xj" ‘ 


Integrating (A. 4) over all space 


a /Vi dv . iti 

dt J 2 J 2 3x. 

nr " 


dv 


(A. 3) 


(A. 4) 


(A. 5) 


If div u "= 0 , 


jl j r u i u i 

dt J 2 


dv 


(A. 6) 


Now look at difference form of (A. 3) 


_6_ Vi _6_ 

5t 2 _u i 6x. u i u j 

i J 


(A. 7) 


Summing over all mesh points 




6 u i u l 


O 1 1 V" 0 

6t 2 “ "I- U i 6x u i u j 

i J 


(A. 8) 
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The RHS of (A. 8) has to be zero as in (A. 6). Therefore du^u^/dx^ must 
be devised such that RHS summation of (A. 8) goes to zero thus conserving 
the total kinetic energy. The difference formula for du^u^/dx^ depends 
on the type of mesh and the number of neighboring points used. 

The following is a fourth order energy conserving scheme using five 
points in one direction. Following usual convention 


f * " ^ [ f + f) + f (x - 

f x = I [ f(x+ f> - -f>] 


(A. 9) 


Note that 


f - f x + 0(A 2 ) 
9x x 


- 1 : V ~ + °“ 2 > 


3f 4 -j-X 1 T-2 x 4v 

■ 3 £ k - 3 £ 2x + 0W > 


(A. 10) 


m (f i-2 - 8£ i-i + 8f i + i - W + 0(a4) 


Now fourth order momentum and energy conserving form of (A.l) and (A. 2) 
are 


3u 


i , 4 ,-xj-X4. 

— + 3 


1 .-2x 1 -2x4. 

Xj ' 3 Cu i ju j j) 2 Xj 
2x. 


3 ( Vx* - 


(A. 11) 
(A. 12) 


To show that these are indeed energy conserving, work with the expanded 
form of (A. 11) . 
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Sa i 

fit 


4 ,-x-x. 
3 ( "i" >* 

-(u y v y ) 
3 U Jl V 'y 

4 ,-z-z. 

3 ( V >2 


l/- 2 x- 2 x. 
3 <U 2 U >22 

l,- 2 y- 2 y. 
3 (u i v > 2 y 

1/-2Z-2Z. 

3 (u 2” >2* 


- 0 


(A. 13) 


where Is either u, v or w . x (A. 13) gives the energy equation. 

Then sum over all mesh points. 


6u„ 


■ 


-£ - T < "f“ 2> ‘ ) 2*} 

+ | V and w - component |1 

[ u ii{t ( V>x + Vx + u< “?x } 

+ . . . J 

• 3 { ( “2">2* + u( “«>2x}] 

+ v, w, - component 


1 ( , . 2 x , - 2 x . — . 2 x 

' 3 ^ <u 2 u> 2 x + u 2 u 2 x + 'W 1 !. 


^ u £ u £ \ f4 -x 1 -2x*| . T 4 
-E— | [I U x ’3 U 2xJ + |_3 


v y 


fl ? . i*v 2s 1 ! 

L 3 z 3 2 z J j 


Iv 2y l 

3 V 2yJ 
(A. 14 ) 
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The first summation on RHS of (A. 14) is zero for periodic boundary 
condition. Therefore (A. 14) becomes 


6'u, 


2>j. it ' -£ 




div u in 4th order; (A. 12) 


(A. 15) 


So if the numerical divergence of velocity is zero, the RHS of (A. 15) 
is zero and (A. 11) is indeed energy conserving. 

The accuracy of the LHS depends on the time advancing scheme and, 
as mentioned earlier, the Adams-Bashforth predictor method introduces 
very weak instability on computational mode. To show that (A. 11) and 
(A. 12) are really fourth order accurate, work with a typical term: 


4 x-x. 1 .-2x-2x. 
T (v u > x - 3 <v u ) 2x 


* IE j (vu) (k+1) * (vu) (k-l> + u (k) (v (k+l) - v (k-l) 

+ v (k) (u (k+l> - "(k-l) 1 | -h | (vu) (k+2) - (vu) (k-2) 

+ U (k) (v (k+2) " V (k-2) } + V (k) (u (k+2) " “(k-2) 5 

Substitute RHS by Taylor expansion as usual. After some algebra, it can 
be shown that 


4.-x-Xv l.-2x-2x» 

3 (v u - 3 (v u >2* 


3 , , A’ 

H (VU) - 60 


f \V , V . V 
(vu) + V + u 


(A. 16) 


Therefore it is indeed a fourth order scheme. An extension to higher 
order differencing scheme can be done on the same basic idea. Overall 
fourth order accuracy is obtained for advective term by using second 
order energy conserving scheme for the Leonard term. 

For convenience, fourth and second order advection term is 
recapitulated below. 
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Fourth 



^ itlK scheme: 

cd et energy c ° nse 

8 i,.)4 

. X 0>,u) + K ( “l 81 

j - 3x 1 , 

+ ^%'3 f °^ 2V V 

U-x-x, .k^2xr 0( 

+ 1 3 CUl “ 1 , (u t ) (V-1) ) 

1 , + u U“l 5 (W« 

e pt l 1 . , r N - ^ u i u) t*“ 2 > 

a I 

. ,\+ “il “Ck+« ’ “°‘' 2) 

+ u((Vo‘ +2 > - (Vd-J 

+ v ( C u i^ (51+1) 

t . r \ N - 0> t v) (t-D v 
A \ U u i v) (ii+i) 1 

1 3 ^ n r v Cu t v) ( i-2) 

\\ -M lUi (2+2) 

* W ai 

v , , - (u t ) (m-1) ) 

+ w U u -s_) (m+1) 

(if \ - C^4 W ) (m-1) 

U (V)(^i) 1 

P L 1 . r N < u i w) Cm-2) 

\\ - -S5 i (u i u1 ^ 

+ l1 ^ , X *,*(»,*» 

( Cu,) (cr f2) ' (*.17) 


+ w 


+ 0 O ) 
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Second order scheme: 


a , N 3 A 3 3 

’K — (u,u.) “ -r— u . u + -5— u . v + u.w 

3 Xj i 3 3 xi 3 y 1 3 z 1 

“ (u*u X ) x + + (-?*), + °CA 2 ) 

■ 44 j'VW) ' (u i u >(k-l) + “ (<V(kH) ' (u lVl>) 

+ U 1 (“(k+1) ' “(k-1)) + (u i v) (i+l) " <U 1 V) (J-1) 

+ v ( (“j) ( 1+1 ) - < u 1 >(l-l)) + u i ( V (H+1) ' v (!-l)) 

+ (u l w) <mH) ' <u i v W» + ” (‘Vtofl) ' '“lVl)) 

+ “i(’ , (» + I) - »0,-l) )j +0 < i2 > «- 18 > 

The subscript is shown whenever it is different from (k,£,m), i.e. 
u (k+l) " u (k+l,Jt,m) 
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APPENDIX II 


FLOW CHARTS AND PROGRAM 


Overall Flow Chart 





INITIATION 

(SUBROUTINES INICON & START) 


CALCULATE EDDY VISCOSITY 
(SUBROUTINE VISCS & VISCV) 


CALCULATE fi ± (EQUATION 3.32) 
(SUBROUTINE VELOC) 


CALCULATE RHS OF POISSON EQUATION 
FOR PRESSURE 

(MAIN, SUBROUTINE DIVGCE) 


PRESSURE FIELD SOLUTION 
(SUBROUTINE PRESS, FFTX, FFTY AND FFTZ) 


ADVANCE ONE TIME STEP BY 
ADAMS -BASHFORTH METHOD 
(MAIN) 


DATA REDUCTION AND PRINT OUT 


NO 


TERMINATE ? 


D 


YES 


WRITE DATA ON TAPE FOR CONTINUATION 
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Flow Chart for INICON 
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nnnnnnnnonnoononf) 


♦ DEC* 

pROC-Pam vain ( T VPtIT.niiTPHT, TAopo, TiPfl'n) 

C****** *«**•*•*• ************************* .♦***♦«**** ****** ************** 

c THIS C0M9P|_S IMF OVERALL SFOHENrE op COmPhT aT I nw, * 

C T^F COHPDT AT TOM IS PpCFfiS^Fo M OSr(_V IN T“F. SmaLL CORE PY TQANSFE&ING * 
c data frp“ Tpp i 4Pf-F cc'°e: »'p«opy, it the e n d of each time step the * 
c SIAIISTKaI pijtKJl TIES 4P£ PPrwTPP nor. at THE FNn OF THg COMPUTA- * 
c T 1 ON vE| r-c 1 TIES aMP THE RIGHT Hano SIDE OF Hp-iEf.TUM FQU A T I ONS ARE * 
C SIPPED E<>r COM JMI'AT TON. * 

I M E C- f «•' TI*'E,2.o,Zf-M,2-“?,2P1 ,7P?.; # 2 ADO I , ZLE SS 1 , ZS A VE , PL ANE 
I . TsTflRT, TEMP 
b Eac x , mlEN.havG 

LARr.E i' M (lA,|E,lM,(;nA,^,tAl l F(|fc,t'>.U)iEHlo,|A.IA)> 

t »1 (1*, Ifc, 16) ,D?( I f-,1*, Im .initf'di, 1 !) ,Pnn9?) 

L ARct , RlJfl 6. t h, 16) 

- , G V (| h , 1 t> . I A> ) , G w f1 » , J <> f 1 6 1 

L APr-t Pl'M*,1E>,1fc1,wT'(16 f 1t>,1M. t »*n*.1b,lM#P)(1f)»l*,1M 
. ,C-I>1{ 1 *>, 1#,, 16) ,RV1 f 1 fs, 1 h, i M ,r.*1 Mb. 1#,, 1 a») 

DIMENSION 1i(Je»lJ>.^i.vn(!*,l6,Sl.’'i(M»,l*«#e)#PU8*l6»5) 

1 ,Cl»H-Y (1 b, I 6,1*) .Sir. fl6, ^,3 ) . F»(ift,16),FI(16,16) 

? , TKFf*,3),TPjr«, 1M.P,rn<S),W«AVE(1(,),NFFTf31 

COhv.On/P a TA | /U, V * * 

COmmOn/pai a?/P 
CONnL’n/DATAJ/SIG#* 

CONhOn /pATAU/O'JMHY 
ruxi-t'N/L.i icj/rjUji.i, !kn,thi 
COWOA /pATA7/FR,Ft 

rnx «Of./r>A tah/k’wat'e, vfft 

CftHfOl. /PA T AQ/ If AT ,.ft 4 y f L MaY,MHALP,NAVG.N| En,nspEc 

COt.F 1st, /(1UII,* DELTA* *2) 

CC’EF2sr*PFLTA*o,S 

roes Uei , /?u , 

COEFS-1 , /(U,*rFLTA**3) 

COFFHsI 2,/DEl H 
C P E F 7 s i , / ( 1 2 , «PEL T A ) 

COf F7 | si ./( 72.*OEL TA**?) 

COE >»*?./ (3, *OT) 

CUEF 9s“I»r>EI.TA*5, /3, 

ALPhAiPT / fUfl.*PELTA) 

PETa«1 ,/fuS t *PEI.T A) 

CPF B ),n Fpw THE F I P S T TTME STEP si,? AFTERWARDS 
l.'THpPf; s ijn*T/M-3,5 

!.' T Ogt,f A N CONVECTIVE OPWVSTPEAH VELOCITY 

-INITIATION 

I *P I T E « 1 
HTHoOGe3«,5 
M0H*1 000./5,0S 
COF I 3 0,5 

CALI I-nIcon f C^E F ] ,rrFFp,cOEF3,cOEFu,CotF5.COEFfe,cnEF7,COEFA,COEF9 

1 ,COpF, , ,rOEF1u,AI.PHA,COF,OELTA,rtT,CONST,i|FR,NHRlTE,GAHMA, 

2 TStaa-t, TE»-P,N«nOEL ) 

NF Il7*1 

kOSsI , / (HFP**? 1 
C.AHKi ,S S r,AMMA**P 
r.C0FSGAT'»*A*CPEF7 
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one* 


Gt> 8 C Af-MA/ 6 ', 

Ri)ac,4f“A/D£LTi 
602*0. 5*G0 

COFr t hB(f>FLT»*.‘JAVG)**2/?tt. 

CUFFUMfoEFl 1 /i:FIT»i** 2 
CDF. F 7 l *C (IF F7 / ( 6 . *OELT A ) 
f»F TAsr’oEFf/u, 

TO t v* 1 . /nulh, 

---rC'NFUT A TTOS ? TAFTS 

00 30 c m«e«TSTART,TEMr> 

10 continue 
c.. ---COMPUTE FOPY VISCOSITY "F". 

2*1 

ZMbI 

70*2 

7P1e3 

C-.-. -FIF$T TRANSFER DATA fBDH LC*< TO SC« FOR thE FLANFS L*L^AX AnO 1 
ST-M.LT** fUC 1 , I , 1 ) .I'M 1 . 1 .L'MT) ,256) 

S T a L I- I * J (V (1 , 1 , 1 ) , VT t 1 . 1 ,t M AY) .256) 

S m A lL I f I'fl » 1 . n ,*’>( 1 . 1 .1 “AY) ,?56) 

SAAtl.IK ( 11 ( 1 , 1 , 2 ) ,'Jf ( 1 , 1,1 ) ,?5 6) 

SMALL Tf- (V ( 1 , 1 ,2) , VMf 1 , 1 , n , 556) 

SHAI LJN (Ml, 1 .2 ),'■** (l.l,t),?56) 

('0 700 L * 1 , LA' A X 
LPIsLyI 

IF (I. '.FO. IVAY) L° 1 * 1 

c transfer vflcctttfs fop t«f plane*i. + i 

SI-AI t. in (U( 1 . 1 ,7PU ,'JMfi , t .LP1 T .256) 

S M a L L 1 m fVM , l ,Z p ! ),VK( t , , ,L p l ),?56) 
fTUl,),?M),-*(i,i,LP) T • ?56 ) 

C THEM “F rAM r r ''-'PUTF FPOV VISCOSITY AT THf p L* n E"L 

r.n TO (?n,25) M'-noFL 
20 CALL VTSCS(ZP,7M,7 p l,7,COEF;>) 
r,o to 30 

?5 CALI VISCV(Z0,7M1 ,7^1 ,7,COEF?) 

30 CONTINUE 

c THF EpPY VISCOSITY ON THE PLA f *E»L is TRANSFePEO to E I** LC* 

SMALLPI.IT (K (1 , 1 ,Z) ,E f 1 , 1 ,L) ,756) 

7SAvE*2M 
7*1sZo 
Z0*7 p l 
?F 1 sZSAVE 
700 CONTINUE 

c. — -rt'MpUTF H(I). THEM STORE t^e * 4 in GU, GV t GW 
7M2sl 
Z *■• 1 * 2 
Z0B3 
2 p 1*4 

Z p ?*5 

SMALL I M run , 1 , U ."*‘(1 ,1 ,LM Y-1 1 ,5121 

SMALL TM (V< 1 , 1, 1) ,V-(1 ,| ,LMAY-n.S12) 

SmalI I* (M 1 , l , 1 ) ,w“f 1 ,1 ,L*AY-n .512) 

S»‘A». 1 1* (P(l , 1 ,3) ,!'►’( 1 ,1 . 1 > ,76l*> 

S*‘ AlL I*: f V( 1 , 1 ,31 , VMf 1 , 1 , 1 ) ,76*) 

.« M A I L I t f fc f 1 , 1 , 3 ) , W :• - ( j , t , 11 , 7 6 ? ) 

S A 1 L I * ' ( * ( 1 , I , 1 ) , E ( 1 , ! , L M * y 1 , ?S 6 ) 

S* Ai L 1 K (*(1,1, 2), F(l, 1,1), 512) 

ZLESS1 «l 
7*2 

;acpi*3 

owg^al page IS 
° F POOR QUALITY 


N 



no 750 p I. A N*E* 1 ,1. 

C*Ll vfLOC (ZP, 7*1 ,Z*'2,?oi ,?P2,?,7LESS» »7APM »BETA,ePEF6,C0EF?. 
1 !)T,COF,GAH'A,r,A**AS,C.cnf ,GD*,GO,nELTA#Ct») 

SMALL CUT fOt'MVYM , I , n ,G('( 1 . 1 .PLANE) ,256) 

St- ALLOt'TLDUM.fv n . 1.21 ,RVf \ , \ , PL »mE 1 , 256 > 
sMALLn»Tfru*"*Y r t , j .3) i . i ,pla*e) ,?56> . 

? y AiLpiiT fSIGft , ! . 1 >.nt M, 1 ,PL**'C) >256) 

smallol^ mr,u,i »2i,n?n ,i ,pu**ej*z^) 

C P|»P«E$SI'RF f OH 1 K J PU T } 0 M h» «e»n BOTTOM 

•sr«A|_LO"T (OLIM-.Y ( 1 , 1 , u ) .run ( j , t ,plA*F.) ,756) 

?M4i LOUT f Ol'MMvr ! . 1 f *) »GV1 ( 1 , l .PLANE) * 256) 

SPALLruT f OI<mmv( I , t ,6) ,G«m () . 1 , PL* *E) >756) 
call oivr.rE (7o,7"i,l>-?.7P' ,?p?.rnEE7) 

SMA L LO'lT (OliMMVfl, ),6),Gf 1, 1, PLANE 1,256) 

7SA vfciZLESSl 
7 L E 5 S 1 * Z 
ZsZaDOI 
7aom*Z8avE 
L*PL* k, ET7 

TF (L ,r,T. L^AX) L*L-L*AX 

FMAJL tM f<( 1 ,1 ,ZA0P1 1 F(1 . \ ,1.) ,256) 

7SAvE*ZM2 

7f2.7«j 

7mj.Zo 

70*7 P 1 

7Pt*7F? 

7P?«7SAVF 
L sP| AnFaX 

IF (L ,GT, Lmakj L*L«LMAy 

sm*lli* M.»n,t,zP?>,yMrt,i 

SMAt L ( V( t , 1 , ZP?> , VHf 1 . \ .1 > , 

* v II* <**U . 1 ,7.P?' \ ,0,256) 

760 CONTINUE 

C«. •••*!>' G 1 1 , (? V , G APE c f1 “ p UTEn Aon DTV 'J IS STORED I* G, 

C fOHpUTF OJVfGL'), THE* GET G 

7*2,1 

7 H 1 1 2 

70b3 

?P1 *« 

7F?,5 

stallin' n»n,i,n,Gi*n,t , lhax«u,5I2) 

S h a l L I tJ f Vf 1 , J , 1 ) ,GVf I , 1 ,|. H A *-! >,512) 

S“alL i a (*< 1 , \ , 1 ) ,C.V‘(L ,1 ,1." AX- 1 1.512) 

S H A L L T F- It'd . t #31, Gun , 1 , t >,T6«) 

SMA L L IN (V( 1, 1 .5) ,PVC1 , 1 , t ) ,76*1 
SHALL IN ( *f f , t . 3) ,C-» (1 , J , J ) ,7681 
no 770 PI. A Mf 8 1 , L M A X 

C ALL OlVGCF(7('.7 y 1,Z M ?,7Pi ,7P?,cnEF7> 

Shall. Oi iT (Oi.'m* Yfl , 1 , J , PL A HE) ,256} 

754vE«7 M ? 

7H?«Z»1 

7H»,7o 
7PszP j 
7 w 1bZP 2 
7P?eZSAVE 
l cPlamtaI 

IF (L .GT. L y A J< ) L*l.-L M AX 

SHALL IN fun , 1 ,7P2 1 .Glif l , I,L1. 256) 

?ma l L (V(j,l,Z p 2)»GVf 1,1,1.1,256) 

SMALL!* f«(l,l,ZP2),Gs« : n,J,L)#256) 

770 CONTINUE 

C-— -O(Gu) IS STOPEO I* PH 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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no 7ec 

OP 7«o J»i , JMA* 
no 7^0 Ib J . IMA* 

PUT. J.l)aCnE***G( l.J,l)*P»MT»J.L> 

7 so r out Rue 

G IS TEMPORARILY STOPEP |N PI 
C..». -COMPUTE M 
7 M 2* 1 
7*1*2 
70 = 1 
7P1*« 

7P?b5 

S'i*Ll I v 01(1,1,1 J.G'JK t,l,L*A*-n,5l?) 

S«*l L (v(i,|,|),ev1{|,l,L«AY-|T,5l?) 

S m «l'l I ‘ f *( 1 . 1 . 1 ) ,r.*i Cl , I ,1* A *-0,512) 

S *s A l L I K (IK 1 , l ,?> ,&Ul f t , 1 , n ,7681 
S M *l l !'• ( V(| , 1 ,5).GV1 ( 1 , 1 ,1 ) .76ST 

smrlr ('(l.i ,J?,r.' Ui , i , 1 1 ,7sa> 

no i 7 7 n PI AN‘P*1 , L* A * 

CAU Mvr.f.f (7 0,2*1 ,2*2, 7»1 ^2,0^7 > 

S**lU'(.-T (OUmmy ( 1 , 1 ,bl,G(l , 1 .PLAN'E) .?5fc) 

ZSAvEaZU? 

7*2*2* ) 

ZM b20 
Z0*7P| 

7P1 sZp? 

?F ?s7S 4 VF 

isPii^e ♦_! 

T E (L .GT, LMAV) LbL-LPAV 
. SMt L LI f OKI O ,2*2), GUI (1 ,1 ,1 
small R (V(l,l ,ZF2>,GVi (I ,1 ,L ),«(>) 

Smalli " f *(1 , 1 ,7R?),G-1 (1 , 1 ,1 >,?5fe> 

1 770 cnOTJM'E 
C-----CO M FUTF ot 

CALL PRESS (COEFJ.COEFt 1 .C0EF7n 
C.--.-*n* STOkF Pi TP Pi F»OM P« A*D S « I r T G FRO* *1 TP G» 
OP lPPO 1*1 # L M A x 

SMALL tv (DUMMY (1 , 1 , 11 .Pi (1 . 1 .1) ,?5fc) 
small ]m (DUMMY ( 1 , 1 ,2) ,P.M( j , l ,I.),?5M 
SMALL CUT (PU U *V( j t , , j i ,ro , i , l>,?56) 

S m A (_ l OUT ( dummy (l,l,?),Pl(l f 1.D.2S6) 

1«Pfi C C“NT I 

c 

PPJVT 905 
00 « 00 L« 1 * 1 <* 

PRTmT 00!J, L 

p M I m T CO?, (E(T,10,U,I»1,!M 
oOO c 0('T T wi/E 
PRINT R 1 7 
no a o i %» b 1 . i f) 
prTmT ct»,J 

PRRT 90?, (F ( I , J, 1 ci , Tal , 1 61 
<401 rOMUHE 
C 

TAIL P«ESS(COEF3,rOEF 1 1 ,roEF7t > 

C ----- s T ORE ",v,w,p at FTvF STFP TnTFpvaL 
T F (A'STORF ,FQ. 5) GO TO 110 
NST('RE«A'STpPF.*l 

r,n to m?o 

110 MTRE bTTHE-1 

write (B) ATlMf.UM, vm,wm,PM^pj 

FND FILE B 
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NSTn»F«l 
. j ?r\ rONilM'JE 

c— pHS*Gii(n-OP/n»(n 

C.««»»S10kF «hSfK) TP «OU*t'. TM?N aoVan CE 1 STEP 
PP POP C* 1 # L^A * 

J P2el-2 

Ip1*t4t 
LP?«l. + ? 

GO TP f7»l ,7*2, 70S, 7*5, 7*5, 7*5, 7*5, 7*5. 7«5. 7*5. 7*5,785. 7*5. 785, 

- 7 *T, 7 M) l 

7*t L^?eLf4X-l 

|. M|*l<*AX 
r.o TP 7*5 
7«2 L K ?sL k, *X 
r.O TO 785 
7*3 | P2*l 

ro TO 795 
7*9 I P 1 e 1 
LP?*2 

7 PS CONTINUE 

PO POP J*1 , JMAX 
JKpsJ-? 

JRlsJ-t 
J p 1 e J ♦ t 
.JP?*j43 

PP to (7*6,7*7,790.790.790,790,790.790.790,790,790,790,790,790, 

- 7**. 7*9) J 
78* jo2eJ“4X-| 

.tf'lsJvix 
r n t i i 7 o/i 

767 

C-C TO 790 
799 jP2s 1 

GO TO 790 
799 TPtal 
JP2*2 

790 fON>TlMIE 

PO 9^0 1 1 1 , I M A t 
IH?aI-2 
T*1eJ-l 
I P 1 a I ♦ 1 
TPPelt? 

GO TO f 791 ,792, 799, 795, 795. 795, 795, 795, 795, 795, 795, 795, 795, 795. 

- 703,799) 7 

791 K?« Jv.ax- 1 
I*laI*4X 
f.P TP 705 

792 l*i2aI."AX 
GO TO 799 

793 I P2e 1 

GO TO 795 
79a TPlst 
IP?«2 

795 COf.TPfF 

sri- f I , n*GH(I, J,l )->C0FF7*(PM(iM ? ,j, L ).8' > *(»M (T ^ 1 ( J,L).PM(tPl, J,L) 
. l»P v (TP?,J,L)) 

!., 2 T *p V(J.J,L)-CPFF7* (PH n,TP2,L)-b'. »(P"(T, J"t , L 1 -PM ( T , JP J ,L) 

- )-PTT,.iP2,L)T 

. BOI'^U, 3T«0<U.v*,U-CPFF7*(pm(I, J , L M 2 )-»'.* (» M C T , J,l*l )-PM(I,J,LPl) 

- )«P*f T,J,LP2)> 

M«(I» J.L)«l» w fI,J,L)*PT*CCOP*omiMf I,n-n,5*Ru(I.J,LM 
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VUtT » J.L 1 »VM(I , J,LltOT*.frnF*«DU>if 1 , 21-0 .5* tf V ( I , J.LJ ) 

W-*«(l * J,LT« fcM ( 1 . J.ll *0T« < COF * POi»m ( 1 , 31-0,5*«w( I , J,L) 1 

eufI.J,U* n O"*fI»n 

pv(j , j,LTs<5P" M n , 2 ) 

e w u»j#n*' ; iiL^n,3) 

HOC CONTINUE 
C-----TL'PnUL£«T E NEE GY 
TKUeO, 

TKVsO, 

T '< a e 0 . 

ru 23o 

t'H't- fL , llaO. 
pI)I.'*'(L.2T«0, 

R0t'"<L* 31sO, 
r-0 230 J*1,.tf’6X 
nn ?3n lst,l”Ax 
T*liaTkl' + ll»( I , J,L) **2 
T« V*Tn V *N“ ( T , J,L)**2 
T«is S Tm<. ) **?. 

w.)U‘;(L# 1 ls^OUMfL/ t HMH( J , J,i_) 
oOUN(t.?>s B r>l;‘HU,2T + w : f 1 . '.LI 
B'Ouvd . 3>sCD<'NfL# . J.ll 

?jn r (.sNTlMig 

ThII s TkII*TPI V 
THV*TkV*TDIV 
T*w s T*W*TDJV 
mSi.imbo , 

VSUM*0’. 

l*SUM»p. 

PC ?*Jp Lal.L^M 

•«***•• #1 4% 

l • » .’ W •• - t • O «. T * • y . • • | U f 1 / 

v5'J*«sYfiu-*i?Pii“(L»2) 

(L » 31 

?ao rONilNUE 

C 

COEel .5 
PWIvT 910 
P 9 I m T 9o« 

T'<Su M aTMUTkV»T , < , T 

PriU.T ®0°*TI v E ,Tkm,TKV,TkI*,TKSI'M 
RTImEsT 1 E * D T 
D J ST *tiF B *R T T M E 
S*ATl0sE*P(0.oi*5<ll*>* ; MST) 

PPJnT 915, PTI' f ,DIST,SP4TrO 
T«Ol Vsl , /TKSU M 
TKU)«TKU*TK0IV 
T“V/t*TKV*T | <OlV 
T*t«isTK*'*TKOTV 

TKt.(r.V.U'Ml/(TKV*T«M) 

PRINT 90‘>,tM'l,T*'Vl,T'"^1,Tm 
I.iT^*Uon*T l*f. *OT ♦IITMOPR 
limial', 0E4O6/TKM 
I'OVsl ,of*0^/Tkv 
"U w-1 , OE + Ofe/T*Ki 

ii04r3,PE*06/TKSli« 

PW T w T 91 1 , UTN , UO" , "OY , UpW , I In 4 
1 • i* a T K 1 1 * * 0 S 

v*«t*.v*wos 
k » a T K U. * 0 S 
P»*'cTt> SHfUwOS 
PPTMT 9Jf» f l'k, VW,Wf,fJN 

r .».» •dissipation TERHSj LFuNABD, SGS, TOTAL. 


97 


ORIGINAL PAGE is 
OF POOR QUALITY 



ni£‘i«o, 

nsr.<$«o. 

00 2«i L»i,Lm4v 
00 ?«) 

00 put in , im«x 
0L€ v «0LF.^'*0l(I ,J,L1 
0SG?aPSCS*0?n , J,u 

2 <11 CONTINUE 

OTPTSOLE“*OSRS 

Ph- 1 n T <nx, r>L F , n$r;S , UTOT ' 

C • • • • • S H E > E 5 S CHECH 
S*3sQ. 

S*2s0. 

f>0 p*0 Ul.INiV 

1*231-2 

TMsI-l 

TPlsUl 

TP?sU? 

T F (I .l-T. 2) C,0 TO 253 
GO tl- 125 1,2521 I 
?5t TNlel*4X 
U -Psl^x-l 
GO TO 

25? T * ? s L 4 X 
GC 10 PSfe 
253 I 1 3 f * 4 X- 1 

IE It .LT. in GO TO 256 
T2 sT*aX-U 1 
no TO (255.25<n 12 
25 ’< IP?sl 

r.o to 25<s 

P'as i*-isl 

1 P?s2 

25* CONTINUE 

no 2*-o J s 1 , .1 u 4 x 

no 280 Lei ,l*it x / 

ounxanm T*p, J,l ) -8 . * (U M (IF 1 . 0# 1 1 1 * » LI 1 •!.•* M p 2 * J »L) 

SKJsSk3e nuo X * * 3 
!4*2sSk2*OUDx**2 

260 CON'TJmUE 
SH3sSt<3*TOI V 
JtKPsf5Kp«T0IV)**t .5 
Skc.sn 3/S*2 

PR I s T 01 2, SK 
r.o TO (270,2*0) NEILT 
270 N F I L 1 * 2 

GO TO 281 
280 N.p Tl T * 1 

261 CONflMUE 
? V I v T 90 8 
PCIvT P 05 

1*0 <iln Lei* mh 4i. f 
pHInT 90«#L 

PPlN.T 902, C'H-n , 1 0,1 ) , lBl,MH4LF) 

TrI' T 902, (V>-n,10,L). Tsl ,NiH41.P) 

p 8 I w T 90?, (HN r J , 1 0,1. ) , ]»1,MH4LF) 
print 902, (Pni-.iCL,n,I*l ,3) 

«10 CONTINUE 

P R I *»T 902, USUH, VStiM,^SUH 
300 CONTINUE 

WRITE (0) TI*:E,Uni,Vw,wk,PM,P!,Rij,RV,Ruj 

F no FILE R 
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320 
90? 
003 
90 a 
905 
908 

908 


909 

910 
91 1 

9 1? 
913 

915 

018 


917 

918 


COM WE 

PUPnAT (X, MEia.fc.xn 

F 09 fAT f///,X,*EP0Y VISCOSITY AT J» 10 *) 

FORMAT (X,*PlAHf. a *,131 

pnc «• a t ( ///, X , *VEL0CI TV AT .IslO » Ut *1 

fok‘«At (x, 1 8 * , sx , o 1 1 * « j / sum* , f i a , 7 , 5 v, 9 HVt*?/ 3 'jMb , Elu, 7 , 

1 <sy. 9 W'***?/S: |k? s , Eia. 7 , 5 x, «k 1 a { V ?-u 2 1 */ { v 2 +n? ) * * . F 1 a’. 7 , RX , 1 H*) 


F » *• A T f y , M H «*«*♦*** , 


, * * * < 


►»*»*«**•»***•«***♦*** 


• +**•»*****■** 




1 


F < 1 1 ! ■ A T (X,1H«, sx, *TI«E s T F P ■ * , Ta,Sx,*u2 a*. Fla, 7, SX, *V<> s*, 

- 5* , =#, e 1 a . 7 , * si |M a* , E i a . 7 , 19x, ih*) 

F CF •' A T (|H|) 

F(iC„A T (y, ifc* l sx,*l!oi/“a*.Fi?.s.3x,*iin2/i!?a*,El?.^,3X,*M02/V28*, 

- F 1 2 ,5 , 3X , »uO?/fc?a« ,F 1 ?,S, xy , *uC ? /H A vG?« * , E 1 2. 5 . 1 a X , 1 N* ) 

FiV’-'AT ( X , 1 8 * , c x , *8* F « *F 8«-» f O't/9 t 1 3 * , E 1 ? . S , 9 1 X , 1 H* 1 

MWAr (X, l'-*,5X,* l '. , ISSIP4Ffn''l LF0»'A90a*,fcl2.5,33, *SGS* * , E 1 ? . S , 

1 )x, ♦ TOT AL s* , F12.S, o 9 x , !*♦) 

t ! ; ••Ur (X, JNi, 8 * , *FFAI. Tt“ga*, Fs’.a. 5** *01 ST A N C E 3 * » F9.U, 

) 4 I V-* , 5 x , * $ T £ A I *4 CaTTOs*. F 7 , 3 , 53X ; tH*) 

cijc-at {X, Ik,, Sy, 1 ?W|, *«?/„,>**? a , £ 1 4 . 7 , "3, 1 2HV**2/^C**2 a , 

1 Mo. 7, ux. |?k,* 0 /k U M? a .E19.7, >jx, a,El«,7, 

2 8 X . l H • ) 

F t> 9 “ A T ( ///, X, *E9pV V ISC n S 1 T V at PLANEal**) 

F"'«‘:Af (y, *Jt*, 13) 

s T 0 p 

f M> 


*DFC* T'lICC'IJ 

S I .'U 9C’I IT I V f I M r 0 fc ( c OF F 1 , e of f f , f rig F 3 , c 9 £f u , r OEf^ , r OF F fc , C OE^ 7 , 
-ri'FFH,C0F eiS ,C!'F.F1 1 ,COgFl«, »L° H A,COP,OELTA,OT#CnfcST,UTM,*l'!RITE, 


2 p A ’ M A , T $ T Awl , T F M O 

c****************************< 

c THIS s Uj. k OUT I Nf tmTIaTES the 
c T I a i PlfLn IS f»E Si 6 0 4 Tg o , Fpc 
C no Tape at Twg goo of t*-f 


PSPC-Pa*. FOR 

CO k, T T mu£ t ion 
previous eiiv a»E 


Start I UR 
P 0 p 9 1_ g M , 
Rf AO I 


pbohleM 

THE OATA 


THg IMI-* 

STOREO * 


C************************* 

1 ' TgCtfi PL A'vF # TS T A 9 T , TFNO 


r ****** *( 


F t A l NO I V, M 2 , W» , »|8I-K , K , WLFu, h A VR 

I. 4»-r;g PMnfc,U,ifci,(;nfc,i^,ifc),Fflfc,16,1fc),gT(t8,18,lfc), 

1 t'l n 8, 1 8, 1 fc) ,r>?M fc, 1 h, 1 fc> ,Fv(fcO) ,EU1 f 80) ,EM (fc«) ,E 02(8*0 
L A F p. F l.-H n ft, 1 b, 1 8 1 , VO ( 1 fc , 1 b , t ft) ,HS n fc, | 6 , 1 M , urn 8, 1 8, 18) 

- , v I n fc , 1 fc , 1 8 ) , ■• I (1 8 , 1 8 , 1 8 1 

I A R R E RU( 1 8. 1 fc, 1 8) , -»V( 1 fc, 1 fc, 1 fc) ,98( Ifc, 1 6. 1M ,Ft ( 1 6, 1 fc, 1 6) 

. , r.U 1 C 1 8 , 1 fc , 1 8 ) , R' ; 1 ( 1 fc , 1 fc , 1 8 1 , Rf 1 ( 1 8 , 1 fc , 1 8 1 

0 i Mg U.STt'.t. fc'- n fc , 1 fc. 1 fcl ,u, 1 8, 1 fc) ,FR ( 1 8 . 1 8) . P'T C 1 ‘ . ! 8) , T9P(B , 31 , 

- TrI ( 8 , 3 ) , Fpo-t ( 1 fc, ?1 , r,9 ( 1 8) , r. I ( 1 6) , U*AVC( 1 », ) , MFFT (3) 
r 0 v *> 0 N / r A T A 5 / fc 9 , G 1 , T P B , T R I 

C. 0 8 1‘ 0 u / o A T A7/F0.F I 
^l)^'^'0^^/PATA«/^.fcAVg,•'lFFT 

r P MM I'.K / o A T A 9 / 1 * 4X, 1 k AX,I..Miy F,V 4 vr,,N| g.V.USPEC 


C»---»f>.ST4F'T*1 STARTING FPOm TTHf STgosO 
ST ARTS2 continued fro** PREVIOUS 4I|ki 

C--»- - r M A y 3 mA y T 1 -FSM VU*'CFjj I fc' y . 0 T ■, g r, T I 0 N 

1UA yivj i 1 ►•I'M *•' F S8* N|pxngp j,\i y-oicECTTOn 

" AXSMAX T HUH MFSH WU-ippp j N 7«0 ) RgCT I 0 : T 

C- — --TSTaRTsST aft I nr time step 
tEnobenoinG TIME STEP 
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c— — DELTA* -FSH SUE 

C-----r>7«T I > F STEP 
C-----rsHnr>7 L C O^S T A V, T 

C Avr,*.">fcLTA ( AVF OiRfMGl /nFLTAf PESH) 

C---- -A MSP* « I* EoijiT f n k (5.16P) 

C--»-»aA M *‘isSTC AIM PATE 
C....»M*<pf'Et a i FOP SMar.C'OI wSf V MODEL 
C--»»» k r’f J pEl*<t FOf* vo^ticITy “FIDEL 

Mao 7rt3,MSTai-T , f k 'A<,.|Mix, L M4x,TSTA9T,Tp^O,MPPOEU 

»£A!i a, PfLTAjf'TjC# ^Avr., AfiT?n,(irM,G4MHA 

f U i-s*' 

kkC S h sl6 

f>*AlF s k .MpS w /? 

k;hP j 8 >'w A if ♦ ) 

► '•'Is'i-FJW.I 

ctSPS3./(5.«A k ISO) 

TE«ps»visO/3. 

PAMSraSTPTnEPP) 

C C 2 1 « 

TFAf*lMfii,*jMA».i_ M AY 

FACsSpoTfTFAC) 

COE FI * 1 , / f 1 an , *oEt. T A**?) 

ro?Fizr.* f 't. :lta*o.s 

rOEF«S| ,/2<l, 

COFFSsl ,✓(<!, *DFLTA**3) 

COP F es 1 2 . /PFL T a 
C OF F 7 s ] . / ( t ? . «OEL T A ) 
n.'EF TpsCOF.F 7*2. 

COEFFs?, /(X. aDT) 
fOFM:^ ai UiSPPoSHSARS 

■* vTF i U-j’. I »F I (?*}•> /N“«Ai.F 

COtF t 1 *COfcF) 0 

C r, FF 1 2*3*. ) 0 1 522o535998 a;>’, 

4U p HA-r)T/fUft. *of L TA) 

C'V $T*COFFIO/ DELTA j 

C 0 M $ T S * f, n N S T ♦ * ? 

Cl'f F l-isCOEFl?*FAC»CONSTS 
C OF f I F>sC0E F l 2*F iC 
P 1 1 aCOFF 1 0 
PJVsPJlA?, 

CALL 5T4CT CC f1 FF3,COEFn .CFlTA) 

or to n ,) non, 5ooo), ^start 
l Cl'Ffl.O 
mCOa'TsI 

2 >»1.25 
YRseGFMfxO) 

? COMlvUE 

C«-*-»E'.'F MY SPECTP'IF F'ATA 

C-- :.1 l T F ° V *. L 1° TO 1.-. The* „S TVTEPValL op To 6.0 

C--»«»Ff tS Thf FA'FkC.v SPECTRUM Frit) Tkf J.SpTPUolC PART. £Vj JS FO» THE 
C AMSOTQ0 P IC PACT, 
no ) '«1 ,2<i,P 

»• 7 s fi ♦ 7 

Ft»0 a, fEAi( v ‘M 
3 COMlM.jF 

on vo* M*i,?a,P 

“7»“A7 

pf. * rii (i , (EM (F-W) »F'M* M . M 7) 

SP3 COFTIMIE 

0 ppiMii (8E10.O) 

00 F> L*1«L M AX 
oo F> Jsl.JF'A* 
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u o 


no 5 i«i,ih*x 
J.L)*0. 

V«?O,J,l)e0. 

*»<] ,~J.L)*0. 

Ml ( !# J,L)*0. 

VI ( ! , J,L)*0, 

*'U!iJiU*0, 

»U(I,J,L )efl. 

*v< l , J,L)aO. 

»<(I, J,l)eO. 

8 rOWTlM'E 

DO ao L«1,VHALF 

I.L*L 

^3«L-1 

F'3$e*3**? 

no 30 .let, mm i 

J IM'ExsJ /ft 

JJa .JJ*«dex* J“AV 

r.;?a.j»NHAl F 
f’2SsM?**2 
r a ?_o isi , 

1 1 f ; CE * *1 /9 

MsI-khalF 

**? 

kSOPbnI s*^?S* f J3S 
T f (NSn« ,11. 0.1) GO TO 20 
WAV^tpPRT (NSPK ) 

L01v*1 ,/WAVN 

*l?s»M Sa* 2S 

ll” i 1 c .C.T, {■ . i ) •vCu' M*tf 

IF (A<*5(H) O, NHAl.f .AMO. A53( f J?) ,F tt . **H *LP 1 NC0 MTb2 
C„.«.-G£T ErnsifcB A^PLITUOP PK TMp 1mtTI»L FTEl.0 *0 PESC g IBED !*•< SEC U,H 
7 »sCO*ST*wa VM 

F' c EG*x*r. i 

GO TO (310,315.315.315.315,315,315) “«EG 
310 »‘**/0,l 

VM S y.o.i*H 
w 1 *«•* j 

FOSE* 

0*VM*in, 

F.OAaEF 1 (-i)-En l (P) 

F amjso«EM < ; 1 WEOA*y*MO. 

GO TC 320 
315 ^s(v.t.)*2, 

Y H «: jr • J ,«0 ,5* M 

***♦10 , 

VlBKf) < 

r Os? A. (•< 1 >.?N(W) 

FMFBGvsE f, {»0 4Fn*vp*2. 

FDAafcM ( V 1 Ji-FA l («) 

?A*!ISPrF.‘ I (•‘')*FCA*V''*?. 

3?0 ftSefNPflGf *K ISO/(COpF i?*v**2) 
ns«?rcT (PS) 

GSAsE.ft* ISft*«IS(i/(C0EFis*x**2) 

OMAsSOS T (OS A ) 

...«-C*AN‘G? t>- a VE ►.■(.«??» VPCTPP Tp .SATT.SFT N'lMpPIeSL DTV F*FF 

...... 1 ,BP. Amt P3 AB? THE ^OOIFIEO "aVE 

GO TP’ (3 30, 3«0) NCOf-'T 
330 COOTIvhf 

* R 0) CPU *M 

ABGg*PX2*Ml 
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B l* <k«?Gn*ST^(*Rf.2n*COEFT2 

AfiGl*PI1**2 

*hG?*PI2»n2 

e?a(P ,*SJ '*'(4Rr.| )«3!Mr4er,2) )*CPFPT2 

APGl *PI 1 4*3 

4fir.?»pt?«vi 

P5S(P ,*STf,(4PGU-S!^r *BG2n*C0EF72 

P 1 Se° 1 * *2 

?2Se'?*«? 

C3S*kj*»? 

P12S*ciS*f?2S 
»S»« B 1 2S* B 3S 
GO to (33?*3u<n fccn*T 
33S CO'TUUE 

o 1 2sSoc T (P 1 ?S ) 

R l?D Iv»l'. / °1? 

R P T V * 1 ,/SOPTfPSQ) 
c ...-GET 4 | u V Ef. TOP 
C p I B ST CHOOSE WAfcOO - PHT 
3*J0 CO^TlMJE. 

Y y spGe M X * ) 

PHlsTV*COEPt? 

CPH^acpS f PHI) 

SPH) *5IM(PMt) 

GO TC ( P, I l) t CO'/T 
« COMlM'E 

4lac-&2*r°HUwi*03*s.-)tv*SPHn*Qi2i'lV 
*?e( 6 J *C e *T*P2*s 5*P>‘ Iv*sphI)*«i?OI V 
*3*- c 1 ?«®OI v*SPHT 
C CALL R4*^O0 M PhI 
xi) 

PHTsY j«COfF ip 
cPHjscnsfPHT) 

SPH]SSIHfPHI) 

Hiscfc2*rPHl*ewR3*Pn|v*sPHj)*Pi?oiV 
R?8f s j *CPH!^P2#P3*POI V*SPHI )*1? Ipol V 
K3s-Cj?#R0lV*SPHt 

r-o to i 2 
n C OA'T I M i£ 

tNnE*s(YY*o,25)*« 

PH!sO.70S3 <, «?*f2*l>JOEX«l) 

» 1 SSI*m(P‘ j 1) 

4 2scOSfPwi 1 
*3*0 « 

YlSCGE^t *1 3 
T »; 0 § X s ( Y 1 ♦0,25)*« 
PHrsO,7flF3902*(?*iMCiEx-i,l 
PlssI,v(Pul ) 

02sCOS(PHl) 

03*o. 

KCOMTel 
12 COVTIVME 

C DETEP'iIn'E * 4fiO B Tv ECU* T I fH; (o‘.b) 

C P A M P C H T up T 4 

Y3*PC’F'v( x 3> 

THE T *sy 3*COEE 1 2 
C AsCfj: f T U ET A ) 

CP*f I n f THPT4) 
n«i 1 1 . JJ.i L )«fjH #C 4*4i 
VP (1 I , J J.LL )«ijN«C4*42 
W»( I ! , JJ,Ll)*0*'*CA*A3 
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wT (1 1, JJ.LUaO^Ca*** 

! F (NJ ,*E. 0) GO TO 20 
vSlGHstVSliD/iS 

VSJGN**HSfR3) /P3 . < ^ 

kffA*:aCk4*CA*RAKISO 

u W 1 1 , J , LL T t i M , LL ) ♦ hq a f>*u.S IGM 

V I C I ! , J J . L L ) * > T( I T . J J . L L W » M * V S I G N 
20 f ON T 1*1 IF 
30 Cf* v T tMiig 
ilO rOMl*UE 

c K.O“ ThF llPOfTB m a L F CP Tw£ K-SRAff has 'AEF ,j OETERmTNED 

C GFT Tmf TRA^SFORmFO VELOCITY AT TH£ COVJiGATF POTnTS 

C--« — CO^JUGA Tf FORM 
C V3* l TO 7, *1 K N?s-T TO 7 

C w 3sL»l - r j 3 * L ►* 

C «N2ajM 

00 ill 1 = 2,9 

I. Msi.t|(»-2*L 

00 u 1 J 3 1 ,1b 

**< JM9WV7 
jhsj* { t e.2«j ) *m 
ro 41 r = i , 16 
*■ s ( T ♦ 1 *5 ) / 1 7 
T“*i*n8-2*n*H 

UM T-, uRfl.J.L) 

vr(i,j,h 
wen*. JM.lTJa ^CI.J.LI 

VI t T", » M .(.H)S-V| u ,.J,L> 

<J1 CONTINUE 

C *3 = 0 , w 1 a 1 TO 7, *23.7 TO 7 

00 i>2 1=2,8 

JM*]M«-2»I 

00 a 2 J* l , 1 b 

M =(J*151/17 

jm s j+(16-2*J)* m 

I F (J >;©, 91 GO TO U? 

ur ( i M , » i m , i )« u»ti,j,n 

VR(IH,Ji,t 1* VR(I,J,n 
k9nM,jv,| kP(t,j,n 
ukjh.jm.i )*-ui(i,j,n 

VI (I*, Jm, J )a-VICI,J,U 
i"!t!",JK, 1 )=-•<! (T,J.I) 

«2 CONTINUE 
C M 1 SN 3=0 

00 u 3 J = 2,* 

JHsj7i«-2*J 

us- r i , j-, i )* iiRf j , j, n 

VR( | ,J>, 1 )• VRfl , J, 1) 

N*( I #.!*■, 1 )• HR(J ,.1,1) 

•j t n , jh, i is-iij u ,.j, n 

VI ( 1 . JH # 1 ) a. v I M , J, 1) 

*1 ( t # JH, t ) *-* f n»Ji 1 ) 

«3 CONTINUE 

C JNVFR’SE TRANSFORM 

C 7 and y transform 

S I Ons. i , 

00 50 1*1 ,l M AX 

SMALL IN fFR(l, 1 ),U»(1 ,1,L),?56) 





o o r» r> 


SHALL! 1 ' (FI <1 , 11,LU(1 , 1 ,Ll,256l 
CALL FFTxfSIGN) 

call ff T v(Sir,^,cc) 

S^ALLniiTf F«( t , i 1 , u*? ( J , 1 ,L) ,?561 
SHALLOT (FI ( I , n»'.'i 1 1 , 1 ,L) ,?S61 
S“4 L L!V (F°(1 , 1 ) , v»n . 1 ,L),?56> 

SH4LLP' (FI (1 » 1 ) « V ! (I ,1 ,L1»?561 
CALL FFTX(SIK'0 

call ff tycsign.cC) 

S"4LLnnT(rp( j , n , v»( | , 1 ,L ) 

SHALL 01 ' I (F t ( t , n.'M ( I , I .1 1 »?561 

SHALL IN (F»(1 . U .«-U(i , i ,L) ,?56) 

SHALL ( FI ( 1 , 1 1 , .it ( t , 1 ,L) '2561 
CALL FFW(SIR^) 

CALL FFTv(SIC.n,CCT 

S'^AlLOHT (F «f ( 1 , 1 ) , *t>{ 1 , » ,L ) » ?S61 

S.HALLOl'TfF I n , 1 ) , «.l n , 1 ,U .256) 

SO CONTINUE 
C Z TpAMSFOP* 

sh4 L li*>' (* M n , i , n.'-^c t . i , i ) .iiooM 
SHALL 1 V (M< 1 , 1 , I ) ,"T r l . 1 . U ,<10661 
CALL FFT?(Sf5»i # MH,M) 

SMALL OUT («M(|,M ) ,;IH ( 1 , 1 , n . <10061 
SmalLoii? ( h( \ , \ , j ) ,UIM , 1 , 1 ) 

Sh.AlLIN (^(t . 1 , 1 1 ,VR( 1 , 1 ,1 T ,'IOQM 

SHALL I ‘I f Md , l , n , vl ( 1 , t . 1 1 .unOfcT 

CALL FFTZCST&n.HM.H) 

SHALL cut (H'lfi , \ , \ ) , Vpf I , 1 , n , uoQfe) 

SHALLOT (Hfi , t , t) , V 1 M , 1 , n , «0<»6) 

small in (HM ( i , i , n ,«m( i . i , n . «n«»6) 
small in fKfi , \ , n . •! n . i , u ,noo*T 
CALL fFT7(SIF.N,HV,m) 

SHAL lout (WMf 1 , 1 , n .-'(I . 1 # M , «09o) 

SMA|.L OUT (M( 1 , t , 1) , M ( t . 1 , 1) ,'JOOVI 
• . • T h F. I u T T I A L FIElO his "EE'' 1 r.F'iEw* TEO. Th£ FOLLOWING IS TO PRIVT 
rUT Tnformatjon On The GENFPaTpo FIELo 
velocities aoe SThkeit T« uR , v o and nr 
---TURrL'LEnt ENERGY CHECK 
TKUsO, 

TKVeO, 

TK«sO, 

TKIsO, 

00 95 l=t , L“A* 

00 95 Jsl , JMAX 
00 05 Isi, ] m A X 
TKUsTxL'TUPfl, J,L1 »*? 

TkV s Tkv+vn(I, J,L)**2 
Tkw c T K'ui'.’p ( j , j, l 1 *♦<? 

9S CONI ImjE 

TOIv- 1 , / u 0 ^ 6 , 

TKUsTK'JtTOIV 
TKVsTk V*TDIV 
TK» = TKt**TDIV 
TKSIIMsTkiuTKVaTKW 
TKIJpcT nii/TK SOU 
TKVBSTK V/TKStIH 

Trwrstkw/tkSUM 
PRINT-* 706 

PPU'T 70O,OT,oflTa,L,k'AvG 
PRIfcT 70?, TK'.i, Tt> V, TK'v, TKSUM 
PRINT 70?, TKUR, TKVR, TKK'R 
PRTnT 706 
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NT I vE»0 
I 1 5 CONTINUE 
PRINT 601 

utotso, 

V T 0 T * 0 . 

M T 0 1 * 0 . 

no i2o 

P° I vT 7JO , l 

USIIM»p’, 

vsu^so. 

N$IImSO . 

00 116 J*1 ,J*AX 
no 1 16 1 = 1 , T*AX 
USU>'Sli$iJ“*"RO , J«l 1 

VSU^svSff+vw ( 1 , J, l 1 

*SU*srfSl‘v*>*»U , J,L > 

116 CONTINUE 

print 7(1 ?, ( UR { I * 1 O , > . ) , T s 1 # 6 ) 

dpInT 70? , ( VS f 1 , 1 ■) , (.1 , 1 * 1 , A 1 

PRUT 7 u ? * !*.«(!, 10,1.), 181,1) 
pRTnT 70?, USUN, Vsiii 1 , v,sum 
liTC.T sut OT+t.'SU* 

VTOTSvtOT+vSU^ 

WTPTSsiTOT ♦ “$0* 

120 CONTINUE 

PRImT 7C?,UT0T , VTOT ,*TnT 
GO TO 130 
1000 COFsl.5 

PSInT 707 
PRINT 706 

»*1»" / 1 1 ,0 1 , UfcL t » ,t 

REAP (10) MI*E,"R,VR,*k,RN,P..I,RV,»W .. 

NTt“EsNTI kl E»l 
PRUT 70S, NT]«E 

PRUT 70h , 

601 FORMAT ( Y, 

PRUT 601 
00 125 L * 1 * L M A * 

PRInT 7 1 0 , L 

PRInT 70? , {IJPTI , 1 0,L1 , IP! , I M *X) 

PRTvl 70?, (VR(I,)0,l),I81,T M AV) 

PSInT 7o?, (wRfI,10,L),I*l,U'AX) 

125 CCNlUuK 
130 CONTINUE 

700 FORMAT (v , * I N I T 1 AL COtvUTTinN, DTs* , F 1 o ,« , * ORl.T as* , F 1 0 , H , 

1 * C**,F7, a, 3*, ‘AVERAGING grips*, F/j, i , * DELTA*) 

701 FC'RnAt fV,*Lc*, 13) 

70? FORMAT (v,«(Et6,7,») ) 

,C3 *"'->•* T (1013) 

70 u FORvAT (OE20.1U1 

705 format (V, *C 0 NTI‘IIIF() a T T I STF° 8*,T<*,/,/) 

706 FURfAT (x, 1 *.******»♦*«♦♦*****“**““****« *********** 



707 FORMAT T1W1) 

710 FORMAT (7,*B|.ANfs*, I3T 

711 FORnAT (*,* INTTTftL CUNfl TT ON* , / , * , *nTs* , E 1 0 , <J , * PELTas*, E10,<l 
- ,♦ <>s* t F 7, a, « I IP**, FJP.a,/) 

Ril TO 60P 0 
5000 CONTINUE 
C0F*1 .0 
©OOP CONTINUE 
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RETURN 

END 


*DFC* PRES* 

SUBBOljtlVE PSESSfciEM.COE*! 1 »CnrF7n 
C************** ************************************************** ******* 

c SOLVE POTSSi)*') EUl'ftTlii'- aY four If* TB i-'JSFuf m s-'ETMon. * 

1 oi (l a, 16. i m #r?n s. i **, i m .pmo-u ,e*'i f j,o).em f * o ) ,pr'i.»‘ , ( i a. » 

2 Pnfl?) 

LARGE IIMY J6, 16, 1M , v*( »*. j *, «r,1 . A. 16, 16) * >) I M h , t * , 1 6 ) 

• , VI M 6, 16. 1 M .M f 1 h. 1 6, 1h> 

LARGE 9U( \ 6, 1 1>, 1 to ,!."/( 1 *, 1 A, 1 6 ) .»*f i *, 1 1> , j A ) , Pi (1 h, \ 6, 1 *>) 

, , Gl.'l (1 b, 1 A, 1 o) . GV . ( i e , t 6. 1 SI ,G<! { 1 A, U, 1 6 ) 

ni^e^sTo*- mmm #*, jh, i >,) , ►,{ i a, i *», i *) 

OI*E n SIOn f»(IMk),F!(tA,tb'#Ti.:5(M),Tci((!,3),RPfl,‘1,8I(lA)' 

, ,MvAVP(1M,MprTf3> 

C0M“On/DA T a 5 /Go ,G l . TRI 

CO^OH/Di T A7/f o, P t 
COBmOn/OAT AH/v*aVp, mFFT 

CO»*HON/0it/i«)/I''A)( . J"4V ,l m * <.mmalp* k, A vG,N|. £* , -*3 PEc 

c PORkABO TRANSFORM 

c FORWARD transform e 4 c u P L»-JF. aftpr T34NSP0RM srn’E PR A FJ TO r, % Pm. 

SlGMB+t . 

00 20 L*1 .Lmax 
SB 4 L I IN ( P # ( I . I \ . C M . I . I 
CALL PPTVfSJG") 

C ALL PPT v f S IG K, ,C r -’tP 
SMALLoi!TfPPM,n.°-'(1 ,1 , L ) » ) 

SBALLOUT (FI(l,l),G(t,1 ,L),2SA) 

20 CO*T I NUE >' 

sba l lin (w-o , \ , n , 0v *n , i , t > ,oooa) 

SMALL IN (Md, J,1 ), 0(1,1 ,1 ),«nP6l 
CALL PPT? (SIC-M.mv.m) 

8BALL0UT »M“f ! # 1 , n ,P‘*(1 ,1 , 1 >,«ORA) 

SmalL CHIT (M(!,1,l) # G(1 # l,j),<ie041 

C-«»-»GE T TpABSPORMEr p R anh pi stored jv r, aNo pm 

NPlC f JHALPfJ 
00 2»0 Lsl.L^AX 
mmsl/ 9 
M S MM*1 6f1 
ARGjecOEPJ 1*(L-W> 

ARG2SAPf,l*2, / 

ARG3 s A°Gl *3. 

arga«arg»*«, 

WAVELsC0S{ARG<*)tl6.yfC n S(A p pj )-COSf APG3) )*A«.*COSf 4RR2)»6S« 

00 200 JsJ.JMAV 
mmsJ/9 
m*mm* t h * t 
ABGJSC0EF|1*(J-M) 

ARG?s a»g 1*2, 

ARG3S 40 G t * 3 , 

ARGaSARGl *Q , 

WAVE J^sCOSt ABG<1)*1 A*. *(Cf>Sf aRGI )-C0S(ARG3) ) ♦*«. *COS ( ARG2 ) -6«? . 

00 20 P Isl , imax 


106 


ORIGINAL paG ^J5 
OF POOR QUALITY 


^GiecOEPiUU-M) 
a»G2*A9G1 *?, 

*PG?s«RGt *3. 

AHGusapg ( *«, 

wAVf leCOSt AHCani fc'.*(COS( apgi )-C<'SCAPG3)T+«jU.*CQS( ApG21«fe5, 

wve(i«Av&TTNAvEJ*i‘-AvEL)«COEF71 

lTES T e!tJ+L 

Tf fA.ssc«v> .LT. D>(inn r,n to 205 
*avf*i./»v 

G(I,J.L ) sG< I , J,U ♦k'AVF 
P"fI»J.UsP* , C!,J,U*«iAVE 
GO TO 20A 
205 G(T,J,U«0, 

P*CI » J,U*0, 

200 CO^TlM'E 

209 CONTINUE 

210 CONTINUE 

C-. — -fxtrapolaTF fop LaSI COpkf.h POTwT 
T*n-pj 
.I sMPl 
LsNPI 

M 1 sNHALF 

K^'2sN».'i -1 

r. f T , J , L )=?.* (G r 1-^1 *G d , M-t .1. )*Gf I , J.^l ) ) 

1 - ( f: ( VMg , j , l l *G ( T , »• <? , L ) *r, ( I , J , N*2 ) ) 

r. ( 1 , J ,1 )s(;n,J,U/T 

p* ( I # ) . Li =2 . * f H •* ( n * 1 , J , l 1 4 f I , »• m 1 , t ) * pm 1 1 , J , nm l n 
1 (\""2 1 4P ‘ ( T ,^m?,lt*p u ( T » J » n >* p ) ) 

P"M # J.L )ePM(I . J.l 1/3 
C tkve&sf t 2 a m s f i. 1 p V i 

Of 100 L s 1 » L M •* * 

S F - * • IT* rFR(i , n . p*(1 . 1,0.256) 

S M ‘LU ( - (F’(M1,G(M,U,?H) 

call m-i vf sir,-' ,c n EF-n 

CALL FFT » (S1GM 

SMALL ru<T(F*M , 1 > 1 #1 ,L) , 256) 

fallout (fjm ,i)i6UiU).?sm 

300 CON'T IMjE 

S H A L L I N (HM(1 , 1 , 1) ,P-( 1 , 1 , 1 ) ,U09«») 

SMALL IP (m(1 , 1 , 11 ,GH , 1 , 1 ),«I096) 
call FFTz (SIGN, *“,.<) 

SMAlLO'.'T (HM( ) , 1 , 1 , 1 , U ,0096) 

SMALLC'iTCHCJ , i , n ,5(1 , 1 . 1 ) ,'4 0PM 

PWImT 0|)ii 

00 9 in (.cl , S 

PFTNT 903.1 . 

PPIvT 901, (P''(1.10,U,Isi,fli 

915 CO^TTLUE 

ooi FURMAt (X.PfEl 0.7,*)) 

903 F(1»vaT (*, *Pl»‘ : F««, 13) 

oou FO»mAt ( 1 HO , *P«E RSi'PE AT .Jain*) 

RETpPN 
E N'O 


*OFCK VELOC 

SI'8POL'Tive VFLOC f 70, 7*1 ,Z‘*2,7Pl Zi ESS 1 . 7 Ano 1 , «ET A , COE. *6. 

1 CrFF7,fJT,COFfGA'"?i,G*MM 4 s,Gc0F.G r *2,Uo*0ELTA,G*') 
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e I***************************************************** ***************** 

c This St'pRr.i.tTT‘:F COMPUTES Hfn In' EGUaTTO*' * 

C***«« ***************************************** ******* ******* *********** 

T'-TFG£<? 7f',2^1,/‘ / 2 < Z p 1»; p ?.7.ZLfSS! ,Z4O0l 

PF*l w , Ni Ew,vavfi 

r* IMf HSU" 1 

- , p (1d,lt,«) ,011**7(1 6.1 k t b) 

rO““nM/CM » l /it, V, w 
CDVMt'.r. /oil ip/P 
f:’“yC«/iJ a TA3/SIG,* 
ro*‘ 1 ' Oi'/oi Ti ti /nitMHy 

rnHwOw/rv 4 1 4 0 / J m a* , JH A X ,1 -4 4 * ( , Mi VQ , N( £W , i»SP£C 

HDtvSl»/'5fc TA 

.ITE$Ts.1‘<ax»1 

ITPsTsT*'**-l 

TO 210 Ort , JHAX 

.tjsj 

.1*1 a .1.1 

J p 1 =J*1 
.1 P?aJ*2 

IP .RT, 2) fin TO 3ft 

r.o TO fjft,20).j 

10 l“2sJ‘AX»l 

.IMs.ImA* 

t;o to ho 
20 JmJ-.Thax 
00 TO bO 

10 TP (JJ ,U. JTf.$T) GO TO bft 
JlSjJ.jHAX+2 
tl* v It i f 5 v'l i -3 i 
U(, J*2s l 

GO TO bO 
bO JMet 
l p ?s? 

bO CONTlM'E 

vs(J-7.5)*0tLTt 

y» s>»r>P. t TA 

VPSV*OFl,T t 

no poo 1 Si , imax 

j v 2sI-? 

r*M el-1 
IPlsUt 
TP?sU? 

IP (I .GT.P) GO TO Oft 
PO TO (7 ft, *01, I 
7ft 1 h? s ]hav.1 
T h | s I m a x 

r n to 

*0 I M ?s I *" A X 
no TO 12ft 

IQ IP (I „l T . JTP-ST) (JO TO 120 
Itsi*TV4V*2 
r.o TO (lOlldlM, II 
ICO IP?sl 

00 TO 1 ?1 
1 10 JPM si 
I p 2s? 

120 rO«Tl»-'UE 

X=f T-7.51*OH». TA 
XT'S y»|>F L T 4 

xPsvtoflta 
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o o r* o 


u component 

.-AOVf r T Ipw CP^POKEWT 

I'M . J. Z?)*u{ I® 1 . J.ZO)-<»n Hi , J,Z0)%Uf ?"1 . J.Zfll 

risrwcct, ipi , j,z<n-i'(Tit ,j.zo)i 

St*$U' ! (T,*»,Zin«(i.f joi # J,ZO)l 

S J *S1 *•'( T , .’“l , 70)*V ( ! ,.JP1 , 7 p > -• I r T , .]n , 70 1 *V( I , f~i , Jo) f 1 , J, ?0 ) * 

-n.'n ».iPt , z <!)•<.' (i, j v i ,zniw "c r. j. /••.»)* rvf I. J p t # z*)-p( r* J"t .23) ) 

S1*S J ♦ '•'{ t #G2Pn*- f 1 , J,7°l '-I'd . t,TPU*«M . J,20)*f 

."(I, J,7P1 J,? w l ) ) flif 1 , .1, 70) ♦{ *( I, I,zn 

Si»*t'(TP£. J. J. 2 *) in )♦*'<!, 

I'.f Id? . 1 . T(\A.n/ti>.3 l.7«\u;in t . >At 7iMailif «9. t . Till \ 


*PVgC :>|n,»S.U2,*S? 

R»*M (]>'?, J.Zn»P.*n>r 1*1 . J,7P).:if !?t 

52 = i>l 1 . J"?, zr > »R,«f i'f 1 , ,i-M , 7ni-n f I , jo i , Z<\ ) 1 -u ( T , ,«P2 , Z<) ) 

53 a 6r0P*(Y*S2-x*S1 T,J,7!>> 

AOV.ECcifVPC. ♦Rr7v*S3 

C--»»- p t'S0L ViPLt SCdtfc (.HRS T|fP“ ft contort T £<?H CftMPfiNfrwT 

S3lsU f TP?, J # 7 o •) * * i ft I, Jnwitf l*>1 f n n»2, J, 20 )-u f T,J,znn 

R J1 sSji ♦ SJ f I°i , J, tP?, J.ZOl-'-'U, J,Zo) ) 

S31 aSji ♦ i * ( IP 1 , I'M ,70) *V ( ?P1 , t p l ,701-IU TPi , J' l , 2 0 I * v ( I ° 1 , J’ 1 1 , Z<1 1 
$3lsSjl *V ( TPI ,,i,7r>). fU(l°I , 1P1 , 7 M-IM I®1 . J'M ,£«)) 

SilsS31 ♦•‘(IPI ,J,Z<n*{''(t»1 ,.|P1 ,7n).Vf ?Pt»J M l,Z«n 
3t *S?1 7P| , J,7P( )*-(1 p l . J, Z p t >-H(T p 1 ,.1,2” 1 ) *'<r T P 1 , J.Z*) ) 

3 s R 3 * * f ’ ^ , ’ * * * ( • ' t I “ ! * j , " ? i )«••*( I P i » I # c ' ! l I ) 

3 1 a S 3 i ♦ ‘ ( j r j # j # 7 i-i > , ( f £ p , , j # 2 p j ) . <■ ( I P ) , ! , ^ j ) ) 

S3?s--if !»•?, r «•?..» , T“i , J , 2 a ) * ( ' ? f i , J , ? o ) - ‘ » f T v ’ ? . J,ZO) ) 

.‘ ! 32sS7?*"( I m i , 1"?, J,Zn)) 

? Um ,.t p 1 ,70) *V(|^t ,.t w l ,‘X0)-'J(TPi I^i , J>"t ,Z0) 

S 5 2 = S 3 ? + v ( J *i i , , 1 , 7 (> > * f ;.l f j m j ,JP1 ,70 )-U(Ip1 . J~M ,10 Ji 
S32 sS 3? ♦Ud- t , J,2/n*<v(T* 1 , .7P1 , 70 >-v ( 1mi ,,f‘M , 2 01 ) 

S32sS??4"C TmI , J, 7p* ) * •>( i“l ,.t, ZPi l-'CT-M » j, l' t (fM . J, Z M 1> 
S32*S3?*f- ( Tut 2rt)* (M( T M J » J. 2Pn*Mf T'M # T# Z?H ) 3 
R3 ?sS3? + ii( t»i , t, Zo)* f *• f i*i , f, 7 pi >-«n , ,r, 2;*m 

Sal sl>( fPt » JP1 , 7«>> *ii( TP) ,.*p) ,7n)."( |mi ,.tP, , 2n1*H(r"l , joi ,70) 

R«i ! aSu 1 tun » JPi , 70)*f "r jpi , yo\ )»iu mi . J P 1 # 3 

SUlsSul ♦ *» c t ».?P1 , ?«>*»•* (jpj , jpj ml . JPJ , 2ft n- 

S«i*saj +'.i<!.JP2.Z<M*v(It.?p£*?fl)+vrt'jPt t Z0)«('!(T»JP?,ZCl- 

,i'( I , J, 7Pn*t'f T,.)Pl,7X)*fVM, JP?,7')>-V(T» ».2o)) 

5«1sS<H ♦iif J, JPl , 2 pi )*••«(!, JPi , 7Pn-"(T* *Pl . Z'\ )*-( I. JPl , 2^1) 

SU1 s S<M *> ( I , JP1 , 70 ) «<<!f j , ,fO» , f T , fPi , Z‘>| ) 1 

s<il sS'j 1 +iif T , JP1 , ?0 ) * f ‘ < r , JP 1 , 7P t )--'(! . .TP 1 , 2*1 ) ) 

Ri|?s"( TPt t Jmi , 7<> Hmi n>! ,.t“i , ?o>.n( tvi , i tf i , 7-1) V»(IrM , >'t , 70) 
i!«?s.Sa ?♦!!(! ,. 1 .- 1 , 2<n* fur r p 1 . 7C)»nf I Ml ,J'<J , 2P) < 

S‘i 2 *SCP 4 UU » J”1 ,zm *( "( 1 PJ ».1''1 , 70 )-•!( l*'l , J M , zn ) ) 

RU2 8 Soi?-ii(T,,TM?,7oi*vf i, .1“?, 7.TUv(T,.J-l ,70i*{UfI, J,Z0)-i'f I, j-?,20) 
1 ) 

R«?s.Ro? 4H( r,.fi , zc ) * ( Vf T, j, 7 P)-Vf I , J*‘2, Zi) ) 

R«?sSa2*. , ( I ,.tv 1 , ?oi )* •» c j ,. 1 w i , 7P1 i-'J{ I , v<\ , 7 M ) * * fl , JM 1 , Z* 1) 
$t>?gSu?*-( I , JM , 70 1 * tur T, .CM , ?PJ )-"f I , T.-M , 7: l ) > 
s«? s SflP*!'i t , v i->j ,7P)*r-;u, j m ,zpn»-in,.i.M ,7 mi ) i 

R51 aJitUPi , J, 7Pn*‘t< JPl , J, 7°1 >-"f T*' 1 , j»7PU*i'U«l , J, 7Pn*<l( l, J, ZPt) 
.* (UfiPU J. 7P1 )-Uf I mi , j, 7 PI )> 

RS 1 eS5 ! i"f I. J.ZPn*ft.'( IPI , J , ? p 1 ) -IK T 1 , .J, ?P ! ) 

• )*i'fI»*IP1 fZPI )*V(f,,)Pl ,7 P 1 )-'!(! ,.r !,7Pl)*^(r#J^l,7' , f) 

S5»sS57 + vf I, J,7P1 i*fiir t,.ipi , 7 pi j, j-aJ , 7° 1 ) ) *i!r j , 2PJI * f v f I # JP t 

, . ?®i )-v(i, j'M ,?pj n 
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$51x551 ♦»(T,.l,7o?)*-*f!,J,ZP2i*«cT,J,7Pn*(»i(I»J,ZP2)- 

,U(T, J.7P))*U(i , .1,7911. fu ( r ,.j,ZP?x--'( r,J#7n) ) 

S5?s"(IPl , J,7*» ) *11(1*1 .J.Z'M i-"(I?M, J,7‘M )*U(lMl , J,Z M 1)+Ud, J,Z*M) 

,*(u ( ipi , j,z*i i.uumi , j.z^i n 

552x552 i"(T,J.Z u n*f'J(lPl,J.7^n-U(!»'l, J»7“l) 

. )♦••( ! . JPI ,Z*tn*v( T . J*1 . ?*M)-"( I, J u l ,7*1 ) *v( j , J«l ,Z'M ) 

$52 r S5?*''U, J, 7*1 )•(•>( T, J»1 . 7*1 . Z-*l) >*'J(I #J.Z M 1 >*(VU. 

,JP1 .7* 1 > - v C 1 * .»- 1 .z^m 

S52sSS? -“(T , J.7'-?1*“U, J. Z'2)*«f I . J.ZM )*(un ,J, ZO) 

»-i> ( I » J . Z' 1 ?) ) ♦ nf I , J , 7 m 11 * f 'v ( T . J . 7 '' ) • l * ( 1 . J . 7 1 ' 2 ) ) 

S6si.*( I ‘’I . J. 70 l * ; iU°l . J. ..)* Z*) *0M“1 , J, ZO) 

Sbe$t*n( T , *(!>( )P) ,.),2f'i-\*f T‘») , J,2Ml 

SfcsSb ♦ (!( t .,t. 70) *(‘J( JPI # f ,Zo)-|i( T«1 , J. ZO) ) 

Sfcxsb*i'fT,jPl,?0)*v(T, IPI . 7o)-n( T, J*1 , 7<‘l * V( I , 70) 

$e.esb*V/( T , J, ?o) • (!lfl . jpi . ZPl-'if T . J-l , ZO ) ) 

Sbssb ♦ J.?M*<V(T,»P1.70)-V (I, JM.ZO)) 

s<js «;*«♦<•{ r,j,zp« j./p i )-"(!. j.Z*n**rT. ‘'rl,.r,zi')*(u(i.j 

Z*i )-n( r ,J. 7* i ))♦"( t ,J. 7ft ) * f~M .J . Z*P - <if I ) > 

R5 Sol. vs -0 ,*» (S? 1 ♦S12+S0 1 ♦ 5'i?*S5l *SS2 -*.*S<>) 

Pe SOL v-uf SPL V**< F "i 

Sir -tP* f v( IP?, .1, 2 0 5 « l. 1 f t , 7 * ) )♦ v * (U( 1* 1 , J®1 » ZO )-M( 1 1»1 . J‘il . 7 rt )) 
$2s-*** fr( I , J, 7« >«■»(?»?, J, 7o) )*v« (n( T*t , .fpj , Z())-ii( r*! , .JM , 70) ) 

S 5s-** f Ilf JPI , joj , 2ft I.i.if VM , IPI , 7ft) IfYP* (ltd , 1P?,701-U( 1 , J, ZO n 
SUx-X* f »>f JM , J‘ l , 7 •> )•!.( T'M . .PM , ?<■ ) )♦ Y**(f( T , J, 70 i.uf I ,.**«?, ZO ) ) 
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S51 rSsi ♦ ,(T . J. 7 p l 1 * fir {P| ,.l, 7 Pn-!i(T*ii , J, 2 Pi) 

.!♦» f I, JPi ,7Pn*vf T, JP1.7P D-'U,.’ 1,7P) ,/P|) 

sSRi + wr t ,.J. 7Pn *f *’f T ,.IP1 . 7°1 i-*n. J‘*i . ZPj > )*«f I , J, 7P| )*(V(T , JPI 

, ,ZP 1 )-v ( I, J ** 1 , ?P 1 ) 1 

S5) s SM ♦•fr.J.2P>l*-rf,.l,ZPZl*-tT,J,79M*(vfi,J,70?)- 

( 1 , J, 7M ) ♦•* ( I , J, 7 e n« c~( I» J,Z e ? '-«(!» J# z«) ) 

^5? s *(iPi , t, z**i >*i*r ip) ,.i, z*^i t^-i ,j#7i , n*"( I'M ,j»7 M n » iji t , j, ?-in 
,*(•**( I P 1 , J i 7 ,J ! 1 )•*•'( I •' I , J , 7 '■* 1 1 1 

SSPeSSP ♦* ( f ,J,ZMl*f IUP 1 ,J,ZM 1 -"(T 1 l,.j#znn 

, 1 ♦«!(!,. IP) , 7 . •*'1 )*V(T ,.1P1 ,2'MI-I* (T ,.I "1 , 7‘ l )*v ( T , J w 1 , 2" 1 1 
SS2sS5P*v( 1 , J,7M ) * (*> ( T, JPI , 7'M l-»f I, J’M ,Z V 1 t , J,7Pl)*(V(I , 

, jpi , z^i )-v ( i , j*» , z-n i 

■S 5 ? a S *i ? -wrr, J,Z M 2) *•"(!» J. Z M ?) ♦*- : ( I . J#7 w n*(i'(I» J.ZO) 

, . ■< r I , J , 7*’ ? 1 1 ♦ •> r I , J, 1 1 * f » (j , J , 7^ ) - 1 ( T , J , 7 % ? ) ) 

S**>‘r I “I , J » 7 •))*•» f IP1 , J, 7P1-- :ri«t , J, 7Pl*ilf I'M » J,7d) 

‘rsexf.'i T, J,7 r ‘l*{> r’P) ,.i,70l»*' p*l ,J,Z.ol ) 

PfesS«> ♦..rT»J,/oi*r i(TPl , i,Z3)-!)(Ti*1 ,J»70U 

Sf>ssr> 4 .'*’f T , JPI , 701*>/( T , jpi ,701- *•'(! , 1”1 #70 )*«(!» »"1 ,70) 

$*z$b+ V( T, J, 7" ) * ( ■'( i , J°1 , 7 0 J - » -f r , J M , Zrt) 1 

?6*S<> **• r T , J « 70 1 * ( \ ( I , IP 1 , ZO 1 «V ( T , J '• 1 , ZO 1 ) 

Sb*S('*i«.'( r , J, 7P1 ) **f f , J, 7P1 l-*-f T . J,ZIM) *-yf T , J» Z *•* 1 1 I, J, ZOlatafl , J 
. . Z p 1 )-w(i , J, 7,11 l ) + i.( 1 , 1, Z01* r \ft , J.Z p 1 1 -•!(!, J»Z'M) ) 
PtSOlvs-o,5*fS31 ♦S32*Sflj ♦S<4?*S c ; l *S52 -r»,*Sfe 1 
Pi SnLvsPFSOuv*i.M.E' J 

Sis n*P« r*(jP2, j,70)-u f j, j,ZP))*v*r *w(Ip1,jpi . ZO )-*«.( ip i.j*<i,zon 
52 b.xm* ( >, t t , i- P, J, 7m 1*v*r , ( t«i , »p) , ZOl-^f p*l , j' l ,Z0) ) 

( jpi , jpj ,2n>-.. ( ,.TP1 ,/0> IfVP* ( ,(l , JP?,70)-«f(T, J»70) ) 

s«s-x*f I P I , JMI ,Z 0 )-:.f ,701 ) ♦ v«** (•„■ f I , J, / O !•»:( I , .}•'■?, ZO )) 

JB1 , J.ZP1 1- •*< 1 *^ \ , J,7P11 1 + YM * f I , JPI , ZB 1 1 I , J *,1 , ZP1) ) 
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St>s-X*rxr I'M, J.Z>*] )-*(TM1 . J,?*1 )>♦*♦<'*( I, JPJ ,Znn — (I,.1“1 , ZmI n 
S7sGCPF*tSt>S?^S3*«!U*S5fS6) 

S9s<-.D;>* ( ( ■* ( I P \ , J , ?P) -V fl -* \ , .1 , 20 ) ) **- ( W ( T . j° 1 # zo ) •* ( \ , J“ \ , 2* J ) * T ) 
PtSADns(S7 + S l n**0IV 
C— — -SI'^-GrTO-SCHE MODEL 

< 3 Ss wUP2,J,70) .w(I,.l,?n)Mi(TPI ,.I,7 p t J.rM > 

S=K{IP),J,Z )*SS 

SSIs ► ( I. J,ZO)-/f l»?,J.7n)+vr Imj, J.ZPJ >-U(T*<) , J,2hj ) 

5*5-* f I *- 1 , J, Z )*SSi 
$<? = COFFfe*S 

ss *(T, jP3,zoi-x(t.j,?0wvcT,.jpi ,?p» j.vn.j'i.rn 

F92N (I,JP 1 ,z )*S 

SS s w(I, J, 70 )-••>( t..J '2, 70 )♦•/(!, , ?P1 )•»/ (I ,J*M ,Z*M ) 

S9as<>m* ( T , JH t , 7)*S$ 

SPSS* 5 *■* ( 1, J.2»Onj )•?.*! •■( 1 ..1, ?*?)•»/ l. J.70)) 

SSsrl'FFS* (So-x ( X , J.ZIESS!)*?'. *( • f J . J. ?o) -“ ( I > », Z M ?) ) ) 

SGSsSS*S<5 

mp>’vY ( r , J ) sOF r a * r ao^fc ARf S^L v* SOStRf 5 Ann) 
r.'jM.AV f t , ..t,*<)sAr!V£cs* i ''fc. , ;AO;'> 

< J r,f i , .1, i )estr , < l , j. n-<>FSPlv*:-f r , J,Z0) 

.« I f, ( I , ?)sSIG f I . J, r , J, 70 i 

?on ro\rI*-!(jf 
?)o rovTlnjt 

Wf TtiRn 
(•Ml 


«. n r » ■< r> * t j«^p- 

*• i ■ c «* • • • * * v<* v r 

sms-ni.T jk f oi^RCf. (/«,?•' i_,z ‘• p, 7 »i « 7 P?,.rOFF 7 ) 

(;*****♦*♦**♦*«*♦♦♦*******•♦**♦**♦*♦***♦♦♦♦**«*♦•***»< 

c TMT 5 SHHhi;»"TT:.«. C At r HI A Tf S TWf r> T VP ur.fc «C F (<I.O,*) 

C T>-f N Th(. v/U'E IS STf'PcD |m Ot.'Mry (f . 

!»*** ******** *************************** * * 1 


FOP 


r * * * i 

rue PL*N€ 


70 * 

* 


1MFOEP 7 0, 7*^1 ,7^?,7 t, l,7 c 2 

fl "SNSli-A "(1 h, 1 6 , , V I j h, 1 6 , 5 ) . M 16 , 1 6 , 5 ) , 0 "**Y f 1 6 , 1 6 , 6 ) 

rO^vOr/D A t A i / u,v,k 
r I)rU'C'A'/l)A7 AO /W'hM» 

AT 49 / [PAX , J‘ A * , L»*AX ,MMA)_F, MA Vr.,MLE'“, k| S p EC 
JTEsTr J>'»X«) 

T TESTsImax-i 
Mn 210 ,1s) , J'lAX 

.IJa.t 

.P»? = J-2 
J'M =J«1 

•' C '.?Sv1 + ? 

IF (J ,G T . 2) r.O TO ?0 
00 TO ( J 0,20), J 
10 ,»M?- JfAX.J 


,1 V 1 S J M A X 
r,n TC> 60 
?n ,!''?= Jma x 
r-o TP 60 

30 IF (JJ .IT. JTEST) GO TO 60 
»1SJJ,JMAX*2 
r.O TO (00,50), .T1 
0 0 vf p 2 c 1 

r.o TP 60 

so JPlsl 
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JP2*2 

60 CONTlM'E 

no ?oo isi , i«ax 
!*2*I-2 
I^tsl-1 
!Plsl< 1 
TP?sI*2 

1 f fl .CT. 21 go in on 
GO TO (70, f'D.I 
7 o 

1 ► 1 * 1 * Ax 
GO TO 120 
hr I M 2 s I M * x 
GO TO 120 

90 IF (I ,L t . 1 7 £ S T ) GO TO 120 
r t s I • T •• 4 v ♦ 2 
GO TO (JOO, J 10), II 
1 f* 0 T P2s 1 

GO TO 120 
1 1 o TP1 si 
TP2*2 

120 C On- T J M.iF 

01 V S "( 1**2, 701 »'-r Tr2,.t, ?nUv ( T , .1 7U)-v ( T , JP2, ?n> 

0 1 V s 0 ! V ♦ > ( I , J , 7 '2 1 • « C I , J , /p ? 1 

Ssit( JP 1 , j, 701-t ( I >1 , .i,7n>*Vf T , JF1 , Z01-'/(T , J?*1 » Z0> 
S*S*-(I,J,2P1)»-(T,J,Z*in 
oi>-s01V+f».*S 
OU1"“V(T, l,F)S r ''TV*CGEF7 
20 0 GO^TlM.'t 
210 rOK'TtM'F 
iv t T c P k 
FNO 


*0FCK' START 

SUPRCijT^E ST 4 ^T TCGEF 3 ,C 0 f - F 1 1 ,gflTA) 

C************»******«**«**********,****«»** 4 ************** 4 ***********a* 

C T*<TS S UP P 00 T T \ £ T*'JTTaT£1 T*f r.r"<$1 am T S FOP fFT POMTpvPS, * 


.*•*♦♦**»*,***> 




n T *• gMS I ov T Rp C o , 3 ) , TP I ( » , 3 ) , r.« n (, ) , r * I ( 1 91 , ( 1 », ) , NFF T ( 3 ) 

rOMMpM/OAT A^/GC f r,| , TRP.TR I 


COPr.C.M/f; 4 r AM /* . A vf , -.’FFT 

CO^mOa/OA T 4 9 / I P 4 * ,..?*’ 4 > ,1 . **a * , ►•H 4 I F , V A VG , "‘I. E V , *SPgir 
DAT 4 t>.*AVF/l ,9, S, 1 3,3,1 1 ,7, 1 G , 2, 10,9, 1 0 . 'J, 12,6,1 bt 


C0FF3=1 , / 1 6 , * * 3 

rot F 1 1=3'. 1 ‘ilSR2^S3SORP/« , 

->r r T l 1 1 SF 

*.FFT(2)=« 

*F FT <3'=2 
00 30 .1=1,3 
TF P ( 1 , Jlst , 

T F I ( 1 , JlsO, 

30 CPvTtk'iif 


no uo r= 2 ,a 

PsF L OaMT )-J . 

H = o«(?f,FF I j 

t f r ( i , i iscosra) 
TK T ( 1 , 1 )s-SIN(P) 
an C O^TlA'i/fc 
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oo so is2,a 
Ln*T( I )-i . 
es2.«p*co&Pii 
TffPU ,21sCPS(P! 
THI (I,2)s-StNfe) 

so CTK'UMiF 

Hbo, *fnFM ) 
TSP(2,^)scr.'S(P) 
TSt (?.. 

PtfiifTs, 

F *0 


*OFC« FFTx 

S:.meO"T j fftx f.s tr,f ) 

C********* A***************************.********************************* 

r FAS* fnnWTfP T«A‘SFO c m 1 f. v »n ( Cc r T I i'm « 

C*»************»***********«*»***»*»*»*********»*******************»**** 

0 1 F Pf 1 1 ►•) »F) f 1 6, 1 kl , T JOf «,3) . T !7 1 f e , 3) , r,s ( 1 M.O T ( U> 

. , * * avf M i>) , » FF 1 ( $) 

cnFur^ /i A i as /(;»?*<" i , i for 
COPr-V'-./ht m/fo,f t 
f O* *»C'-/P /. T A *■ / >>• A v f , • -F F 7 

CC**-On/da r AO/l H a* , .!< AX ,L |>’A * , i-'-.A, F ,:,'AVG,H E\V'SP£c 

TF fStr.v .LI. r.,) or. TO 3 
rn ? j= 1 . J"AK 
oo i ] 8i,i* 

FJfT.Jlso. 

1 CO ,J T 1 M'F 

? COf.TlMlE 
3 CO*T 1 AiME 

00 IPO Jrl.JNK / 

JPej 

PO 20 mm « l , 3 
TE*r. = 0 

INCRSNFFT (MM) 

T F so 

S CO^TlMlg 
TSTaKT= 1*IP 
TEPO*lST API *1 V .CR»1 
vcl 

CO 10 IsJSTART , JEHO 
JPs J+IAC& 

goov t sff n , .n *fp n p , tp ^ 

0l)"«2sFi n , J)*r t (JP, JP) 

0 * 3sF" f ] ( .M »f » ( ! D » .1 P ) 

ROOusM ( 1 ,J)-F I ( TO, JP) 

GlUif-Sr r,(- 1 1.-.3 * t pf ( v . i ..-niiv <( * r w j ( v,, -*f' ) *§ 1 on 

finir*!Aisf;oi , '-3*rpT (‘ 

F»(I . 

F) ( I , Jlsr.OO ''2 

FRf lP,.’PAs(5Ml*’t 
FI (IP.JF )*r,OUM*> 

VSM+l 

10 CO*!TJW|»F 

IF (IP .1 T, I M A x 1 do in 5 
20 COK-ilAUF 

C F I F T m 1,2,3, FTC. 

60 00 70 1st , IMAA .2 
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IPsI^l 

GOU*lsFftf 1. J)*PRC!P,.IP) 
GOUt?sFTri»J)*FIfIP f JP) 
GOl'vJsFW f 1 , J)-FC ( IP, JPJ 
r.ci*^«sn f i, jj-fiiip, JP) 
r.Rfj lsru-cri 
r, i u iscru*? 
r.Ki p ):crii''j 
G I U p )sGM.J- m u 
70 C<“MTfMif 

C t.E t t 4 s nsf uwi’j f value OR^EPEri 

no » n T = i,Jf!A« 

I*** «4Vf( I J 
FR< t •* # , 1 \ * 6 R < I T 

p/i CO»'T f MiF 

100 COMIm.E 

BEltfFSvi 

FND 


*OEC h FFTy 

StfMPCi.lT 1 V F. FPTV 

C*********************************************************************** 

C PAST pOhPjPw T& A v5F •>?■.• j.-j v.fifSfoTTO'i * 

C****************** **♦*****♦***♦♦«*« *♦»**•***♦*♦*******♦♦»*****♦********* 
IX k 5TC- p^t. * *■ . 1 *>1 1 » f }*•• 1 a \ . i »p t #. 3 ) . To j ( H , 31 , OR f 1 h) ,r,I ( 1 b) 

. ,L.,iVf flM/FFTHI 
C<» n.-.'/i.iA raS/r-P, o r , T»G, TPI 
C(.u*f<U>"/n»TA7/F«,F i 
CU“f P‘ /l'-4 T( F ✓'.!> iV(-' , . ;CFT 

C 0 H li(V,/p 4 TaV/ 1 VflK, ji, 4 » >L v A )f # N !t , iL e,M ft VG,fJ ( £M, Kf SPEC . 

c V-TRAMSP n«M 

nr ?00 1 a) * JFA* 

IPs I 

nr i 2 o F u s 1 , 3 

lEKtssn 

Jfvf&BVKFT ( ) 

.TPeo 

105 CONTINUE 

TSTaPTsI *JP 
TEMn=ISTA«TtT*r«-l 
F> s 1 

no tic Jsistart, I£vd 
.JPsJU TNC» 

gouk 1 *r«f i» JUfM jp,.iPi 
r,DI'M 2 sP I ( J , J)*F T f IP,.JP> 
pot’ k JsP«n,.j)»pp{ip,.?Pi 
nout-fsP i n , jt-f t < tp,,»pi 

r.M'M5sr,ut>?»3*-TWM f *, m'- > .r.n'f 'O *trt >*SIgn 

Gl.>i.»'*0sGC-t." 1 3* Tt T (O’, M”i *Sj r. v ;*r,nti«u*TRH ff> ,'iMi 
PR(! , JlsGUUtj 
PI ( J . J)sr,nUM? 

FFf IP, »P)sGCUN5 
PI <l p , JPUGOfft-F 
««»»♦ 1 

uo r ontImif 

fp (vtp ,IT. JOAX) GO in i o*> 

120 CONTINUE 
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X 



nO 1 7 0 JaJ , J"AX,? 

JP*J*1 

GOU*l*FP(I. JT*FP( tP# JP1 
r.oi’ fc ?sp i n , j ) ♦ f i r i » » jp i 
r,oi)H3*rPf r , J)-F*-(|P, 
PDOKU-ej n #.n»FK i p » .IP T 
If (SIGN .1 T , 0.) G 1 ' Tn Ufl 
GP( J)eGDi't'\ 

Gj ( J ) sGC-IimP 
G»< jP)sr,nut '3 
G! ( J P )sGIH |H u 

go to i7n 
ton continue 

gh( j)rr.r>"N» *cr<FP3 
r,H jt*r,ru''?*ctFF'3 

G«(vt p tsGM.'“3*crFF3 

GI(J p )sG'H ,v ‘p*CCEF3 

170 CONTINUE 

c r,ET r-KoecfO SFT 
no i t , r jsi *J M tT 
AVE f J ) 

F«( I. Jf'JsGPC Jt 
FI(I»J v )aGlCJ) 

IPO CONTINUE 
poo CONTTNijE 
RETuKn 
ENO 


•DECK EFTZ 

SUP. POUT I N'E FPT 7 ( S TGU , HM, N ) 


C**««** 
C FAST 


FOiifir JFK 


»*#* + *♦ 

T* A*S Pi. >r ^- 


*++«+«* i 

T ** 7»0J*FfTTnv 




*♦*♦**♦*****♦« #* 


C I M E n s I C *■' M h, 10. 1 Ol , Mf if,, i o, i 0) ,HOnFM »5, 2) ,TPP ( P , 31 » TPT (P, 3 ) , 

-GM 1 0) , l-U 1 bl ,n, «vf M 6T , *FF T (3) 

CC'«m)N/I.'&T4S/GP .gt,to«,t»i 
rC' K HON/OATtF/»KAVF,«.FFT 

CONKC-v/r.iA Ti 9/1 Ni * , AX ,1. “4 T, NHAl. F, NA VG, N| F N.^SPf C 
00 ?00 JSl.IPAX 
no 20r JS1.J0AX 
no i n u «i«3 

TEN'OSO 

TNC9=UFFT(PM) 

IPSO 

1 C 5 CONTINUE 

IST4»T=1 alp 

TENdbiSTAPTaTKcP-J 

*»t 


no 110 Is I ST apt » I EUO 

l.Psi* I Nf-P 

OU^l *H M (I , J,L ) *HP( I , .1,1 P) 
nUMpauf 1, J,l)4»-(I,.I,I.P) 

PUHJ*h«( T,.!,L)»NN(T,.I,LP) 
nu*«l«*( I , I , 

rH)N 5 = ni.'> 3 *TpP(p» >, n)»piima,*t 9 t fN ,vm)*stGN 
nuMo c r>i)M3*Tm )*sir.N*ou»T 0 *TwP( M ,“M) 

F*(I , J.LIsOUKl 
M(I, J,L)cCUF? 
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H"(I# J,lP)oOUM5 
^CI , J,lPlsOU*«> 

M«M+1 

1 10 CONTINUE 

IF (LP ,LT. L*AX) GO TO 105 
120 CONTINUE 

DO 1 70 l*l»LMAX,2 
L P *L ♦ 1 

DU*) sm(I , J,L) ♦ *{ T . J. LP) 

RI'MJ*h(I,J,L)-«{I,J.LP) 

r , t ,.r,L p ) 

hdiivcl. n«f>u*i 

HOor (L,2)s0tjM? 

MCUH (L°i 1 )*PUP3 
MPU M (lP,2)«DU"<i 
1 70 CiJf.TlM.iE 
C GET OPDfPEO SET 

r.O 1*0 Ist.LMAX 
(,l*tK.»4VF (1 ) 

M( I , J,L fc ?SHr-l.'"fL# 1 1 
Mf"( I , J,l * J sf-n'i«' (L , 2) 

1«n f cm T lwi.it 
?0« CONTINUE 
BE ThRn 
E MU 


DECK vjSCvOct 

5 i.l R P l< U T 1 n E visrv ( fr., 7 1 ,7P1 , 7 , Cf'F p 2 ) 


C****o*****«****“***‘***** 

c THIS SUPPm.'rm CAt-CLiTrs T Mf 
>***•*•*»*»«******«•*»+* 

T^Tp GE m 70*2 f 1 , 7 P I .7 
PEAL * . MLE f- » f A v G 

DIMENSION <J('<,,16,^)#v/(U,l(!,,57,M(l^,l4,,S)»SlGUi!,,l«i,35,t<(1f!#)6»3) 
r OHi'Of /DM A 1 /U, V , - 



r.t'f'»-'J‘“/DATA3/STG,i< 

CC'M^nv/OaTAP / lf-i* , Jr* A* ,LH«X,».-HA|.F,NAVG,N| t^.MSPEC . 
c rOF.F? = r*0FLT A*0 ,5 
DO 210 .1 Js 1 , Jf. A X 
JHlsJJ-1 
JPISJJ+1 

T F (JJ ,f(l. n jPia.Jf-AX 

tF fJJ , f O i Jha*) ,»P1«M 
00 ?O0 !»1 , IMAX 
TfUM 
IP1 « 1 ♦ 1 

TF (I '.Eft. 15 ThI.tmax 
IF (I . F.^ , !►'*») T P 1 a 1 

SI * (*fl , J p l , 70 )-.. M , J**i , 70)-vf( T , J.J,7P1 )*v( J, ))♦*? 

S?» <U( I , JJ,7 p 1 )-»'f T . .1.1 .7^1 3 -S* ( TP1 , J.T, 70 > ♦ - c I r 1 ,JJ, 70)) **2 
93s ( vr I p l , JJ,70)-V (1M ,JJ,7o 1-0(1 , JP1 ,20-1 ♦ !'( I , J M 1 , ZO) 7 **2 
S<Jss 1 aS?*S 3 

n ( I , >J J * Z 1 aCOEF ? *Soe i ( sa ) 

200 COMTIM'E 
210 CONTINUE 
B£TllHN 
ENP 


119 

ORIGINAL PAGE IS 
OF POOR QUALITY 



<*OFC* vJSf-Si'AG 

«bo«C(.<TJ‘ £ VTSCS(Z0,7’M ,7P1 .?,Cn£P2) 


c TMtS Sbp»n».TTM| C*LC'H_*Tf.5 T*F ENiv >/ 1 SC 'i-S I I v Fy S •• ifino I MS* y ''’DEL, * 



TMfGr? 7 fl . Z 1 ” . Z e l , Z 
WfcH X, "l ?N,> Ayf- 

f f.»r*M0x:/p4 T t J /ti, V , X 

r«.*»*“CM/n4T*x/s rr-,x 

T 4 :# / T l ’i X , .t *■ 1 X , I. « S * , MM A | F , M4 V/G ,f«| £(w , O' S P P T 

C CCFF?=r*rtl.TA*o.5 

.ITF. 5 T 5 .IK 4 x-j 
T t V s T s T •' a » - 1 
DO ?1P J.’sl , Jv.4» 

.MPsJJ -2 
,JM aJJ-1 
,!Pl=.Jj4t 
.JP? = JJ4? 

TP (JJ ,01, ?i Tii in 
r.n to (h,?o),jj 
jo .)► 2sJ> ! **-1 
.» f. 1 sJf-4* 

GO TO *0 

A . «,«•».. tii I V 

e o ■ 

r. o to t,n 

30 IF fjj .LT. .tTEStl GO TO hf< 

.1 x *? 

r.n to (uo,50), ji 
op Jp?«l 
r.n to ap 
50 ,lP2s? 

JPlsI 

to CONTlM'E 

ro 2 op isi,iMt* 

T 1*2=1-? 

1 "• I = 1 - 1 
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